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We  have  developed  computational  methods  for  the  identification  and  optimal  control  of 
distributed  parameter  systems.  Our  work  has  consisted  of  a  theoretical,  an  experimental,  and  a 
numerical  component.  Using  functional  analytic  techniques  we  demonstrate  and  study  the 
convergence  properties  of  our  schemes.  Extensive  numerical  studies  are  then  carried  out  in  order  to 
fully  assess  their  feasibility,  performance,  and  limitations. 

Central  to  our  general  approach  is  the  notion  of  approximating  infinite  dimensional 
optimization  problems  by  sequences  of  finite  dimensional  ones.  Typically  this  involves  the 
approximation  of  infinite  dimensional  state  dynamics  (in  the  form  of  distributed  parameter  systems 
such  as  partial  differential  equations  or  hereditary  systems)  by  sequences  of  finite  dimensional 
dynamical  systems  (such  as  ordinary  differential  or  difference  equations).  In  the  case  of  the 
parameter  estimation  problems,  when  it  is  functional  (i.e.  spatially  or  temporally  varying) 
parameters  which  are  to  be  identified,  the  infinite  dimensional  admissible  parameter  spaces  must 
also  be  discretized. 

In  general,  we  have  used  polynomial  and  Hermite  spline  function,  as  well  as  modal  function 
based  finite  element  methods  to  accomplish  these  tasks.  We  note  that  since  we  are  not  simply 
solving  the  so-called  forward  problem,  i.e.  the  straight  forward  integration  of  the  underlying 
dynamical  equations,  but  rather  are  seeking  to  solve  problems  whose  solutions  depend  in  a  highly 
nonlinear  fashion  on  the  system's  infinite  dimensional  state  transition,  input  /  output  or  parameter 
space  /  output  maps,  the  theoretical  components  of  our  investigations  become  especially  important. 
Indeed,  methods  which  have  been  known  to  do  a  good  job  integrating  differential  equations  often 
perform  less  than  satisfactorily  when  coupled  with  a  scheme  to  solve  a  parameter  estimation  or 
control  problem. 

On  the  other  hand,  we  have  found  our  numerical  studies  to  be  extremely  useful  in  allowing  us 
to  observe  the  limitations  and  short  comings  of  our  methods  and  to  identify  important  directions  for 
future  research. 

Below,  we  summarize  our  results  and  our  progress  in  as  yet  incomplete  but  on-going  projects. 
In  the  body  of  the  report,  we  simply  provide  a  brief  outline  and  broad  summary  of  our  findings. 
The  actual  results  are  discussed  in  detail  in  the  research  papers  which  have  been  provided  to 
AFOSR,  and  a  sampling  of  which  have  been  included  in  the  appendix  to  this  report. 
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1.  Control 


a.  We  have  developed  an  abstract  approximation  framework  for  the  discrete-time  linear 
quadratic  Gaussian  (LQG)  control,  estimator,  and  compensator  problems  for  systems  whose  state 
dynamics  are  described  by  linear  semigroups  of  operators  on  infinite  dimensional  Hilbert  spaces. 
The  computational  schemes  included  in  the  framework  yield  finite  dimensional  approximations  to 
the  optimal  feedback  control  laws,  estimator  gains  and  LQG  compensators.  A  convergence  theory 
has  been  established  and  numerical  studies  involving  parabolic  (heat  /  diffusion),  hereditary,  and 
elastic  systems  have  been  carried  out.  Initially,  our  theory  applied  only  to  control  systems  whose 
continuous-time  input  and  output  operators  were  bounded.  However,  we  have  been  able  to  extend 
these  results  to  apply  to  discrete-time  systems  whose  underlying  continuous-time  formulations 
involve  unbounded  input  and/or  output  maps.  An  unbounded  input  operator  is  one  that  maps  the 
control  space  into  a  space  larger  than  the  standard  state  space  in  which  the  problem  is  usually 
formulated.  An  unbounded  output  operator  has  domain  smaller  than  the  standard  state  space.  We 
have  been  able  to  successfully  apply  our  abstract  theory  to  distributed  parameter  systems  with 
boundary  control,  to  hereditary  systems  with  control  delays,  to  boundary  control  systems  with 
control  delays,  heat/diffusion  equations  with  pointwise  measurement  of  temperature,  beam 
equations  with  pointwise  measurement  of  strain  or  acceleration  and  distributed  systems  involving 
output  delays.  The  computational  schemes  that  we  implemented  and  tested  were  either  polynomial 
spline,  Hermite  spline,  or  modal  function  based  and  were  able  to  handle  reasonably  complex 
problems  when  run  on  a  micro-computer  (an  IBM  PC  -  AT).  TTiis  was  joint  work  with  Professor 
J.S.  Gibson  of  the  Department  of  Mechanical,  Aerospace  and  Nuclear  Engineering  at  the 
University  of  California,  Los  Angeles. 

b.  We  have  developed  an  a-shift  technique  which  can  be  used  in  conjunction  with  schemes  for 
the  optimal  LQ  stablization  of  hereditary  systems.  This  leads  to  control  laws  which  yield  a 
prescribed  degree  of  stability,  i.e.  all  system  poles  to  the  left  of  the  line  Re  z  =  a  in  the  complex 
plane.  Both  the  continuous  and  discrete  time  cases  were  considered.  This  was  also  joint  work 
with  Professor  Gibson. 

c.  We  have  started  to  investigate  and  develop  a  finite  dimensional  approximation  theory  for  the 
design  of  optimal  fixed  finite  order  compensators  for  distributed  parameter  systems.  The  approach 
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we  are  taking  is  based  upon  and  uses  the  Hyland-Bemstein  optimal  projection  equations;  a  set  of 
necessary  conditions  for  optimality  which,  in  an  infinite  dimensional  setdng,  take  the  form  of  a 
coupled  system  of  operator  Riccati  and  Lyapunov  equations.  We  replace  the  infinite  dimensional 
system  operators  (i.  e.  A,  B,  and  C)  by  finite  element  approximations.  The  resulting  finite 
dimensional  system  of  coupled  matrix  Riccati  and  Lyapunov  equations  are  then  solved  using 
effective  and  efficient  finite  dimensional  optimal  projection  algorithms  and  software  developed  by 
Hyland  and  Bernstein.  At  present,  results  from  numerical  studies  carried  out  on  examples 
involving  heat  and  beam  equations  and  hereditary  systems  are  promising.  Further  computational 
studies  along  with  theoretical  analyses  (i.  e.  convergence  arguments,  etc.)  are  currently  underway 
and  continuing.  This  work  is  joint  with  Dr.  D.  S.  Bernstein  of  the  Harris  Corporation  in 
Melbourne,  Florida. 

2.  Parameter  Identification 

a.  We  have  developed  computational  methods  for  the  estimation  of  spatially  varying  material 
parameters  (specifically  flexural  stiffness  and  Voigt-Kelvin  viscoelastic  damping  coefficients)  in 
Euler-BemouUi  models  for  the  vibration  of  flexible  beams  with  and  without  tip  appendages.  Our 
schemes  involve  spline-based  finite  element  discretizations  of  the  second  order  in  time,  fourth  order 
in  space  system  of  partial  differential  or  hybrid  system  of  ordinary  and  parital  differential  equations 
and  the  function  space  admissible  parameter  set  A  convergence  theory  was  established  and 
extensive  numerical  studies  using  simulation  data  was  carried  out  on  both  conventional  (serial)  and 
vector  computers.  The  schemes  performed  satisfactorily.  This  was  joint  work  with  Professor 
H.T.  Banks  of  the  Division  of  Applied  Mathematics,  Brown  University  and  Dr.  J.M.  Crowley  of 
the  United  States  Air  Force  Academy. 

b.  We  have  tested  our  general  approach  for  the  estimation  of  unknown  parameters  in  models 
for  the  vibration  of  flexible  structures  on  actual  expermental  data  taken  from  the  RPL  experiment. 
The  RPL  structure  consists  of  four  flexible  beams  cantilevered  to  a  rotating  hub.  The  structure  was 
designed  and  built  (and  currenty  resides)  at  the  Charles  Stark  Draper  Laboratory  (CSDL)  in 
Cambrigdge,  MA  with  support  from  the  Air  Force  Rocket  Propulsion  Laboratory  (RPL)  (now  the 
Air  Force  Astronautical  Laboratory  (AFAL))  at  Edwards  Air  Force  Base  in  California.  Using 
accelerometer  data  we  were  able  to  successfully  identify  parameters  in  a  distributed  parameter 


model  for  the  structure.  We  hope  to  continue  to  test  our  theories  and  methods  on  this  structure  in 
the  future.  This  was  joint  work  with  Professor  Banks,  Mr.  S.S.  Gates  of  CSDL  and  a  former 
student,  Ms.  Y.  Wang  who  is  currently  a  research  associate  in  the  Division  of  Applied 
Mathematics  at  Brown  University. 

c.  We  have  initiated  a  collaboration  with  Dr.  Alok  Das  and  his  associates  at  the  Air  Force 
Astronautical  Laboratory  (AFAL)  at  Edwards  Air  Force  Base  in  California  for  the  purpose  of 
collecting  experimental  data  for  use  in  the  testing  of  our  theory  and  computational  schemes  for  the 
identification  of  distributed  parameter  systems.  Specifically,  we  have  made  two  visits  to  AFAL  and 
taken  data  from  an  experimental  flexible  5  ’  x  5  ’  alluminum  grid  which  has  been  constructed  by 
Dr.  Das  and  his  group.  We  are  planning  to  develop  appropriate  theory  and  computational  methods 
which  can  be  used  to  fit  a  two  dimensional  distributed  parameter  model  to  the  structure.  In 
addition,  we  are  also  currently  planning  a  series  of  flexible  structure  experiments  to  be  carried  out 
in  the  spring  of  1988  in  the  30ft  thermal  vaccuum  chamber  at  AFAL.  The  purpose  of  these 
experiments  is  the  collection  of  data  for  use  in  the  study  of  thermal  effects  on  internal  damping 
mechanisms  of  composite  materals.  Our  general  approach  will  involve  the  identification  of 
appropriate  distributed  thermoelastic  or  thermoviscoelastic  models.  Appropriate  models,  theory 
and  computational  schemes  are  being  developed  as  the  planning  of  the  experiments  and  the 
preparation  of  the  vaccuum  chamber  and  experimental  structure  continues.  The  primary  motivation 
for  this  investigation  is  the  solar  heating  of  orbiting  large  flexible  spacecraft.  This  effort  is  joint 
with  Dr.  H.  T.  Banks  of  the  Division  of  Applied  Mathematics  at  Brown  University  and  Dr.  D.  J. 
Inman  of  the  Department  of  Mechanical  and  Aerospace  Engineering  at  The  State  University  of  New 
York  at  Buffalo. 

d.  We  have  developed  an  abstract  approximation  framework  for  the  identification  of  parameters 
in  nonlinear  distributed  parameter  systems.  Using  the  theory  of  monotone  operators  and  nonlinear 
evolution  systems,  we  establish  convergence  results  for  Galerkin  finite  element  methods  for  inverse 
problems  involving  broad  classes  of  autonomous  and  nonautonomous  nonlinear  partial  differential 
equations.  This  new  nonlinear  theory  completely  subsumes  the  existing  linear  theory  and  serves  to 
generalize  many  of  our  earlier  results.  In  addition,  it  can  be  applied  to  parameter  estimation 
problems  for  a  frequently  cited  model  for  nonlinear  heat  conduction.  In  addition  to  the  theoretical 


I 


results,  we  have  carried  out  some  preliminary  numerical  studies  on  a  vector  machine  with  the  aid  of 
a  grant  (of  computer  time)  from  the  San  Diego  Supercomputer  Center.  This  work  is  joint  with  Dr. 
H.  T.  Banks  of  the  Division  of  Applied  Mathematics  at  Brown  University  and  Dr.  S.  Reich  of  the 
Department  of  Mathematics  at  the  University  of  Southern  California. 


3.  Publications  Carrying  AFOSR  Grant  Number  AFOSR-84-0393 


1 .  A  Numerical  Scheme  for  the  Identification  of  Hybrid  Systems  Describing  the  Vibration  of 
Flexible  Beams  with  Tip  Bodies,  Journal  of  Math  Analysis  and  Applications,  1 16  (1986), 
262-288. 

2.  Spline-Based  Rayleigh-Ritz  Methods  for  the  Approximation  of  the  Natural  Modes  of 
Vibration  for  Flexible  Beams  with  Tip  Bodies,  (Quarterly  of  Appl.  Math.,  Volume  XLIV 
(1986)  169-  185. 

3.  Approximation  Methods  for  Inverse  Problems  Involving  theVibration  of  Beams  with  Tip 
Bodies,  Proceedings,  23rd  IEEE  Conference  on  Decision  and  Control,  Las  Vegas,  Nevada, 
December,  1984. 

4.  A  Galerkin  Method  for  the  Estimation  of  Parameters  in  Hybrid  Systems  Governing  the 
Vibration  of  Flexible  Beams  with  Tip  Bodies,  (with  H.T.  Banks),  ICASE  Report  No.  85-7, 
Institute  for  Computer  Applications  in  Science  and  Engineering,  NASA  Langley  Research 
Center,  Hampton,  VA,  February,  1985. 

5.  Approximation  Methods  for  the  Solution  of  Inverse  Problems  in  Lake  and  Sea  Sediment 
Core  Analysis,  (with  H.T.  Banks),  Proceedings,  24th  IEEE  Conference  on  Decision  and 
Control,  Ft.  Lauderdale,  Florida,  December,  1985. 

6.  Numerical  Schemes  for  the  Estimation  of  Functional  Parameters  in  Distributed  Models  for 
Mixing  Mechanisms  in  Lake  and  Sea  Sediment  Cores,  (with  H.T.  Banks),  Inverse 
Problems,  2(1987),  1-23. 

7.  Numerical  Approximation  for  the  Infinite-Dimensional  Discrete-Time  Optimal 
Linear-Quadratic  Regulator  Problem,  (withJ.S.  Gibson),  SIAM  J.  Control  and 
Optimization,  26(1988),  to  appear. 

8.  Shifting  the  Closed-Loop  Spectrum  in  the  Optimal  Linear  Quadratic  Regulator  Problem  for 
Hereditary  Systems,  (with  J.S.  Gibson),  IEEE  Transactions  on  Automatic  Control, 
AC-32(1987),  831-836. 

9.  Estimation  of  Stiffness  and  Damping  in  Cantilevered  Euler-Bemoulli  Beams  with  Tip 
Bodies,  (with  H.T.  Banks  and  C.  Wang),  Proceedings,  Founh  IFAC  Symposium  on 
Control  of  Distributed  Parameter  Systems,  Los  Angeles,  CA,  June,  1986. 

10.  Computational  Methods  for  the  Identification  of  Spatially  Varying  Stiffness  and  Damping  in 
Beams  (with  H.T.  Banks),  Control  -  Theory  and  Advanced  Technology,  2(1987),  1-32. 

1 1 .  Methods  for  the  Identification  of  Material  Parameters  in  Distributed  Models  for  Flexible 
Structures,  (with  H.T.  Banks  and  J.M.  Crowley),  Mathematica  Aplicada  e  Computacional, 

2  (1986),  139-168. 

12.  The  Identification  of  a  Distributed  Parameter  Model  for  a  Flexible  Structure,  (with  H.T. 
Banks,  S.S.  Gates  and  Y.  Wang),  SIAM  J.  Control  and  Optimization,  to  appear. 

13.  Computational  Methods  for  Optimal  Linear  -  Quadratic  Compensators  for  Infinite 
Dimensional  Discrete-Time  Systems,  (with  J.S.  Gibson),  Proceedings  of  International 
Conference  on  Control  and  Identification  of  Distributed  Systems,  Springer- Verlag  Lecture 
Notes  in  Control  and  Information  Sciences,  to  appear. 

14.  Inverse  Problems  in  the  Modeling  of  Vibrations  of  Flexible  Beams,  (with  H.  T.  Banks  and 
R.  K.  Powers),  Proceedings  of  the  International  Conference  on  Control  and  Identification  of 
Distributed  Parameter  Systems,  Springer-Verlag  Lecture  Notes  in  Control  and  Information 
Sciences,  to  appear. 


15.  Approximation  of  Discrete-Time  LQG  Compensators  for  Distributed  Systems  with  Boundai  ^ 
Input  and  Unbounded  Measurement,  (with  J.  S.  Gibson),  Automatica,  to  appear. 


16. 


17. 


18. 


19. 


Approximation  in  Discrete-Time  Boundary  Control  of  Flexible  Structures,  (with  J.  S. 
Gibson),  Proceedings  of  the  26^  IEEE  Conference  on  Decision  and  Control,  Dec.  9-11, 
1987,  I^s  Angeles,  CA,  to  appear. 


Computational  Methods  for  the  Solution  of  Infinite  Dimensional  Discrete-Time  Regulator 
Problems  with  Unbounded  Input  (with  M.  A.  Lie)  Proceedings  of  IMACS/IFAC 
International  Symposium  on  Modeling  and  Simulation  of  Distributed  Parameter  Systems, 
Oct.  6-9,  1987,  Hiroshima,  Japan. 


Approximation  of  Discrete-Time  LQR  Problems  for  Boundary  Control  Systems  with  Control 
Delays,  Proceedings  of  IFIP  Conference  on  Optimal  Control  of  Systems  Goverened  by 
Patial  Differential  Equations,  July  6-9,  1987,  Santiago  de  Compostela,  Spain, 
Springer-Verlag,  to  appear. 


An  Approximation  Framework  for  the  Identification  of  Nonlinear  Distributed  Parameter 
Systems,  (with  H.  T.  Banks  and  S.  Reich),  in  preparation. 
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4.  Meetings  Attended,  Talks  Given,  and  Papers  Presented 


Invited  Participant.  Workshop  on  Control  Systems  Governed  by  Partial  Differential 
Equations  with  Application  to  Large  Flexible  Structures,  The  Pennsylvania  State  University, 
Clearwater,  FL,  March  4-  8,  1985. 

Invited  Speaker.  Applied  Mathematics  Seminar,  Department  of  Matliematics,  Harvey  Mudd 
College,  Claremont,  CA.,  January  31,  1986. 

Invited  Speaker.  Control  Systems  Seminar,  Departments  of  Mathematics  and  Electrical 
Engineering,  The  Institute  of  Technology,  University  of  Minnesota,  Minneapolis,  MN,  June 
5,  1986. 

Invited  Speaker.  Conference  on  Control  and  Indentification  of  Distributed  Systems,  The 
Institute  of  Mathematics  of  the  University  of  Graz,  Vorau,  Austria,  July  6-  12,  1986. 

Invited  Speaker  and  Session  Chairman.  Meeting  of  the  Society  for  Engineering  Science, 

State  University  of  New  York  at  Buffalo,  Buffalo,  NY,  August  25-27,  1986. 

Invited  Participant.  Second  Workshop  on  the  Control  of  Systems 
Governed  by  Partial  Differential  Equations  sponsored  by  AFOSR, 

NSF  and  the  University  of  Montreal,  Val  David,  Quebec,  Canada,  October  5  -9,  1986. 

Invited  Speaker.  Control  Systems  Seminar,  Department  of  Electrical  and  Computer 
Engineering,  University  of  California,  Santa  Barbara,  Santa  Barbara,  CA,  October  27, 

1986. 

Paper  Presented.  1984  IEEE  Conference  on  Decision  and  Control,  Las  Vegas,  Nevada, 
December,  1984. 

Paper  Presented.  1985  SIAM  Fall  Meeting,  Arizona  State  University,  Tempe,  Arizona, 
October,  1985. 

Paper  Presented.  1985  IEEE  Conference  on  Decision  and  Control,  Ft.  Lauderdale,  Florida, 
December,  1985. 

Invited  Speaker.  IFIP  Conference  on  Optimal  Control  of  Systems  Governed  by  Partial 
Differential  Equations,  July,  6-9,  1987,  Santiago  de  Compostela,  Spain. 

Speaker  and  Invited  Session  Chairman.  IMACS/IFAC  International  Symposium  on 
Modeling  and  Simulation  of  Distributed  Parameter  Systems,  Oct.  6-9,  1987,  Hiroshima, 
Japan. 

Attendee.  ICIAM  '87,  First  International  Conference  on  Industrial  and  Applied  Mathematics, 
June  29  -  July  3,  1987,  Paris,  France. 


law 


5.  Students  Supported 

1.  Ms.  Y.  Wang,  MSEE,  University  of  Southern  California,  1984,  MS  Applied  Mathematics, 
University  of  Southern  California,  1986.  Carried  out  computations  for  identification  of  RPL 
structure.  Thesis:  An  Inverse  Problem  for  a  Flexible  Structure.  Supported:  1  June,  1986  - 
31  July,  1986. 

2.  Mr.  M.  Lie,  BSEE,  University  of  Southern  California,  Carried  out  computations  for  optimal 
discrete-time  LQG  compensators  for  infinite  dimensional  systems.  Supported:  1  June,  1986 
-  31,  May,  1987 

3.  Mr.  P.  Feehan,  MSEE,  University  of  Southern  California,  Carried  out  computations  for 
preliminary  studies  on  the  identification  of  material  parameters  in  distributed  parameter  models 
for  flexible  structures  using  modal  or  spectral  data.  Supported:  1  June,  1986  -  31  August, 
1986. 

4.  Mr.  C.  Lo,  MSCE,  University  of  Southern  California,  BSEE,  George  Washington 
University,  carried  out  supercomputer  calculations  for  studies  on  the  identification  of 
nonlinear  distributed  parameter  systems.  Attended  San  Diego  Supercomputer  Center  Summer 
Institute,  Summer,  1987.  Supported  1,  June,  1987  -  Present. 

5.  Mr.  C.  Mao,  Sc.  D.  Mathematics,  Wuhan  University,  carried  out  preliminary  theoretical 
studies  on  thermomechanical  models  in  flexible  structures.  Supported  1,  June,  1987  -  31, 
August,  1987. 


6.  Equipment  Purchased 

1.  IBM  PC  AT  and  peripherals.  Used  to  carry  out  many  of  the  computations  reported  on  above. 

2.  AST  Premier/286  and  peripherals.  Used  by  P.  I.  and  students  to  carry  out  computations 
reported  on  above. 


A  Numerical  Scheme  for  the  Identification  of  Hybrid  Systems  Describing  the  Vibration  of 
Flexible  Beams  with  Tip  Bodies,  I.  G.  Rosen. 

A  cubic  spline-based  Galerkin-like  method  is  developed  for  the 
identification  of  a  class  of  hybrid  systems  which  decribe  the  transverse 
vibration  of  flexible  beams  with  attached  tip  bodies.  The  identification 
problem  is  formulated  as  a  least-squares  fit  to  data  subject  to  the  system 
dynamics  given  by  a  coupled  system  of  ordinary  and  partial  differential 
equations  recast  as  an  abstract  evolution  equation  (AEE)  in  an  appropriate 
infinite-dimensional  Hilbert  space.  Projecting  the  AEE  into  spline-based 
subspaces  leads  naturally  to  a  sequence  of  approximating  finite-dimensional 
identification  problems.  The  solutions  to  these  problems  are  shown  to 
exist,  are  relatively  easily  computed,  and  are  shown  to,  in  some  sense, 
converge  to  solutions  to  the  original  identification  problem.  Numerical 
results  for  a  variety  of  examples  are  discussed. 

Spline-Based  Rayleigh-Ritz  Methods  for  the  Approximation  of  the  Natural  Modes  of 
Vibration  for  Flexible  Beams  with  Tip  Bodies,  I.  G.  Rosen. 

Rayleigh-Ritz  methods  for  the  approximation  of  the  natural  modes  for  a 
class  of  vibration  problems  involving  flexible  beams  with  tip  bodies  using 
subspaces  of  piecewise  polynomial  spline  functions  are  devloped.  An 
abstract  operator- theoretic  formulation  of  the  eigenvalue  problem  is  derived 
and  spectral  properties  investigated.  The  existing  theory  for  spline-based 
Rayleigh-Ritz  methods  applied  to  elliptic  differential  operators  and  the 
approximation  properties  of  interpolatory  splines  are  used  to  argue 
convergence  and  establish  rates  of  convergence.  An  example  and  numerical 
results  are  discussed. 

Approximation  Methods  for  Inverse  Problems  Involving  the  Vibration  of  Beams  with 
Tip  Bodies,  I.  G.  Rosen. 

In  this  short  paper  we  briefly  outline  two  cubic  spline  based 
approximation  schemes  for  the  solution  of  inverse  problems  involving  the 
vibration  of  flexible  beams  with  attached  tip  bodies.  The  identification 
problem  is  formulated  as  the  least  squares  fit  to  data  of  a  hybrid  system  of 
coupled  partial  and  ordinary  differential  equations  describing  the  dynamics 
of  the  beam  and  tip  bodies.  The  resulting  optimization  problem  is  infinite 
dimensional  and  as  such,  necessitates  the  use  of  some  form  of 
approximation.  The  schemes  we  have  developed  are  based  upon  the 
construction  of  a  sequence  of  approximating  identification  problems  in 
which  the  underlying  constraining  state  equations  are  semi-discrete  finite 
dimensional  approximations  to  the  infinite  dimensional  distributed  system 
which  governs  the  original  identification  problem.  Our  study  includes  both 
theoretical  convergence  results  and  numerical  testing. 

A  Galerkin  Method  for  the  Estimation  of  Parameters  in  Hybrid  Systems  Governing  the 
Vibration  of  Flexible  Beams  with  Tip  Bodies,  H  Thomas  Banks  and  I.  G.  Rosen."^ 

In  this  report  we  develop  an  approximation  scheme  for  the  identification 
of  hybrid  systems  describing  the  transverse  vibrations  of  flexible  beams 
with  attached  tip  bodies.  In  particular,  problems  involving  the  estimation  of 
functional  parameters  (spatially  varying  stiffness  and/or  linear  mass 
density,  temporally  and/or  spatially  varying  loads,  etc.)  are  considered.  The 
identification  problem  is  formulated  as  a  least  squares  fit  to  data  subject  to 
the  coupled  system  of  partial  and  ordinary  differential  equations  describing 


the  transverse  displacement  of  the  beam  and  the  motion  of  the  tip  hoodies 
respectively.  A  cubic  spline-based  Galerkin  method  applied  to  the  state 
equations  in  weak  form  and  the  discretization  of  the  admissible  parameter 
space  yield  a  sequence  of  approximting  finite  dimensional  identification 
problems.  We  demonstrate  that  each  of  the  approximating  problems  admits 
a  solution  and  that  from  the  resulting  sequence  of  optimal  solutions  a 
convergent  subsequence  can  be  extracted,  the  limit  of  which  is  a  solution  to 
the  original  identification  problem.  The  approximating  identification 
problems  can  be  solved  using  standard  techniques  and  readily  available 
software.  Numerical  results  for  a  variety  of  examples  are  provided. 

Numerical  Schemes  for  the  Estimation  of  Functional  Parameters  in  Distributed 
Models  for  Mixing  Mechanisms  in  Lake  and  Sea  Sediment  Cores,  H.  T.  Banks,  and 
I.  G.  Rosen. 

We  consider  distributed  parameter  models  for  vertical  mixing  in  lake 
and  sea  sediment  cores.  Finite  dimensional  approximation  schemes  are 
developed  for  the  solution  of  associated  inverse  problems.  The  schemes 
permit  one  to  estimate  temporally  and  spatially  varying  functional 
parameters  which  appear  in  the  parabolic  partial  differential  equations  and 
boundary  conditions  constituting  the  models.  Theoretical  convergence 
results  are  established.  Numerical  findings  are  presented  which 
demonstrate  the  potential  of  the  methods.  An  example  involving  the 
identification  of  a  depth-dependent  mixing  parameter  based  upon  volcanic 
ash  data  from  the  North  Atlantic  is  included. 

Numerical  Approximation  for  the  Infinite-Dimensional  Discrete-Time  Optimal 
Linear-Quadratic  Regulator  Problem,  J.  S.  Gibson,  and  I.  G.  Rosen. 

An  abstract  approximation  framework  is  developed  for  the  finite  and 
infinite  horizon  discrete-time  linear-quadratic  regulator  problems  for 
systems  whose  state  dynamics  are  described  by  a  linear  semigroup  of 
operators  on  an  infinite-dimensional  Hilbert  space.  The  schemes  included 
in  the  framework  yield  finite-dimensional  approximations  to  the  linear  state 
feedback  gains  which  determine  the  optimal  control  law.  Convergence 
agruments  are  given.  Examples  involving  hereditary  and  parabolic  systems 
and  the  vibration  of  a  flexible  beam  are  consider^.  Spline-based  finite 
element  schemes  for  these  classes  of  problems,  together  with  numerical 
results,  are  presented  and  discussed. 


Shifting  the  Closed-Loop  Spectrum  in  the  Optimal  Linear  Quadratic  Regulator 
Problem  for  Hereditary  Systems,  J.  S.  Gibson  and  I.  G.  Rosen. 

In  the  optimal  linear  quadratic  regulator  problem  for  finite  dimensional 
systems,  the  method  known  as  an  a-shift  can  be  used  to  produce  a 
closed-loop  system  whose  spectrum  lies  to  the  left  of  some  specified  vertical 
line;  that  is,  a  closed-loop  system  with  a  prescribed  degree  of  stability. 
This  paper  treats  the  extension  of  the  a-shift  to  hereditary  systems.  As  in 
finite  dimensions,  the  shift  can  be  accomplihed  by  adding  a  times  the 
identity  to  the  open-loop  semigroup  generator  and  then  solving  an  optimal 
regulator  problem.  However,  this  approach  does  not  work  with  a  new 
approximation  scheme  for  hereditary  control  problems  recently  developed 
by  Kappel  and  Salamon.  Since  this  scheme  is  among  the  best  to  date  for  the 
numerical  solution  of  the  linear  regulator  problem  for  hereditary  systems,  an 
alternative  method  for  shifting  the  closed-loop  spectrum  is  needed.  An 
a-shift  technique  that  can  be  used  with  the  Kappel-Salamon  approximation 
scheme  is  developed.  Both  the  continuous-time  and  discrete-time  problems 
are  considered.  A  numerical  example  which  demonstrates  the  feasibility  of 


the  methcxl  is  included. 


Estimation  of  Stiffness  and  Damping  in  Cantilevered  Eulcr-Bemoulli  Beams  with 
Tip  Bodies,  H.  T.  Banks  and  I.  G.  Rosen. 

We  develop  finite  dimensional  approximation  schemes  for  the 
identification  of  spatially  varying  material  parameters,  i.  e.  flexural  stiffness 
and  viscous  damping  coefficients  in  hybrid  models  for  flexible  beams  with 
tip  bodies.  Our  schemes  are  derived  via  an  application  of  spline-based 
Galerkin  techniques  to  the  conservative  form  state  space  representation  for 
the  coupled  system  of  ordinary  and  partial  differential  equations  and 
boundary  conditions  which  describe  the  dynamics  of  the  system.  A 
convergence  theory  is  briefly  outlined  and  a  discussion  of  our  findings 
based  upon  extensive  numerical  studies  carried  out  on  both  conventional 
and  vector  processors  is  included. 

Computational  Methods  for  the  Identification  of  Spatially  Varying  Stiffness  and 
Damping  in  Beams,  H.  T.  Banks  and  I.  G.  Rosen. 


A  numerical  approximation  scheme  for  the  estimation  of  functional 
parameters  in  Euler-Bemoulli  models  for  the  transverse  vibration  of  flexible 
beams  with  beams  with  tip  bodies  is  developed.  The  method  permits  the 
identification  of  spatially  varying  flexural  stiffness  and  Voigt-Kelvin 
viscoelastic  damping  coefficients  which  appear  in  the  hybrid  system  of 
ordinary  and  partial  differential  equations  and  boundary  conditions 
describing  the  dynamics  of  such  structures.  An  inverse  problem  is 
formulated  as  a  least  squares  fit  to  data  subject  to  constraints  in  the  form  of  a 
vector  system  of  abstract  first  order  evolution  equations.  Spline-based  finite 
element  approximations  are  used  to  finite  dimensionalize  the  problem. 
Theoretical  convergence  results  are  given  and  numerical  studies  carried  out 
on  both  conventional  (serial)  and  vector  computers  are  discussed. 

10.  Methods  for  the  Identification  of  Material  Parameters  in  Distributed  Models  for 
Flexible  Structures,  H.  T.  Banks,  J.  M.  Crowley  and  I.  G.  Rosen. 

In  this  paper  we  present  theoretical  and  numerical  results  for  inverse 
problems  involving  estimation  of  spatially  varying  parameters  such  as 
stiffness  and  damping  in  distributed  models  for  elastic  structures  such  as 
Euler-Bemoulli  beams.  An  outline  of  algorithms  we  have  used  and  a 
summary  of  our  computational  experiences  are  presented. 

11.  The  Identification  of  a  Distributed  Parameter  Model  for  a  Flexible  Structure,  H.  T. 
Banks,  S.  S.  Gates,  I.  G.  Rosen,  and  Y.  Wang. 

We  develop  a  computational  method  for  the  estimation  of  parameters  in 
a  distributed  model  for  a  flexible  structure.  The  structure  we  consider  (part 
of  the  "RPL  experiment")  consists  of  a  cantilevered  beam  with  a  thruster 
and  linear  accelerometer  at  the  free  end.  The  thruster  is  fed  by  a  pressurized 
hose  whose  horizontal  motion  effects  the  transverse  vibration  of  the  beam. 

We  use  the  Euler-Bemoulli  theory  to  model  the  vibration  of  the  beam  and 
treat  the  hose-thmster  assembly  as  a  lumped  or  point  mass-dashpot-spring 
system  at  the  tip.  Using  measurements  of  linear  accleration  at  the  tip,  we 
estimate  the  hose  parameters  (mass,  stiffness,  damping)  and  a  Voigt-Kelvin 
viscoelastic  structural  damping  parameter  for  the  beam  using  a  least  squares 
fit  to  the  data. 


We  consider  spline  based  approximations  to  the  hybrid  (coupled 
ordinary  and  partial  differentia]  equations)  system;  theoretical  convergence 
results  and  numerical  studies  with  both  simulation  and  actual  experimental 
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data  obtained  from  the  structure  are  presented  and  discussed. 

12.  Computational  Methods  for  Optimal  Linear-Quadratic  Compensators  for  Infinite 
Dimensional  Discrete-Time  Systems,  J.  S.  Gibson  and  I.  G.  Rosen. 

An  abstract  approximation  theory  and  computational  methods  are 
developed  for  the  determination  of  optimal  linear-quadratic  feedback 
controls,  observes  and  compensators  for  infinite  dimensional  discrete-time 
systems.  Particular  attention  is  paid  to  systems  whose  open-loop  dynamics 
are  described  by  semigroup  of  operators  on  Hilbert  spaces.  The  approach 
taken  is  based  upon  the  finite  dimensional  approximation  of  the  infinite 
dimensional  operator  Riccati  equations  which  characterize  the  optimal 
feedback  control  and  observer  gains.  Theoretical  convergence  results  are 
presented  and  discussed.  Numerical  results  for  an  example  involving  a  heat 
equation  with  boundary  control  are  presented  and  used  to  demonstrate  the 
feasibility  of  our  methods. 

13.  Inverse  Problems  in  the  Modeling  of  Vibrations  of  Flexible  Beams,  H.  T.  Banks, 

R.  K.  Powers  and  I.  G.  Rosen. 

The  formulation  and  solution  of  inverse  problems  for  the  estimation  of 
parameters  which  describe  damping  and  other  dynamic  properties  in 
distributed  models  for  the  vibration  of  flexible  structures  is  considered. 
Motivated  by  a  slewing  beam  experiment,  the  identification  of  a  nonlinear 
velocity  dependent  term  which  models  air  drag  damping  in  the 
Euler-Bemoulli  equadon  is  investigated.  Galarkin  techniques  are  used  to 
generate  finite  dimensional  approximations.  Convergence  esdmates  and 
numerical  results  are  given.  The  modeling  of,  and  related  inverse  problems 
for  the  dynamics  of  a  high  pressure  hose  line  feeding  a  gas  thruster  actuator 
at  the  tip  of  a  cantilevered  beam  are  then  considered.  Approximation  and 
convergence  are  discussed  and  numerical  results  involving  experimental 
data  are  presented. 

14.  Approximation  of  Discrete-Time  LQG  Compensators  for  Distributed  Systems  with 
Boundary  Input  and  Unbounded  Measurement,  J.  S.  Gibson  and  I G.  Rosen. 

We  consider  the  approximation  of  optimal  discrete-time  linear  quadratic 
Gaussian  (LQG)  compensators  for  distributed  parameter  control  systems 
with  boundary  input  and  unbounded  measurement.  Our  approach  applies  to 
a  wide  range  of  problems  that  can  be  formulated  in  a  state  space  on  which 
both  the  discrete-time  input  and  output  operators  are  continuous. 
Approximating  compensators  are  obtained  via  application  of  the  LQG  theory 
and  associated  approximation  results  for  infinite  dimensional  discrete-time 
control  system  with  bounded  input  and  output.  Numerical  results  for  spline 
and  modal  based  approximation  schemes  used  to  compute  optimal 
compensators  for  a  one  dimensional  heat  equation  with  either  Neumann  or 
Dirichlet  boundary  control  and  pointwise  measurement  of  temperature  are 
presented  and  discussed. 

15.  Approximation  in  Discrete-Time  Boundary  Control  of  Flexible  Structures,  J.  S. 
Gibson  and  I.  G.  Rosen. 

This  paper  treats  discrete-time  LQG  optimal  control  of  flexible 
structures  with  boundary  control  and  what  normally  are  unbounded 
measurement  operators. 

The  application  of  recently  developed  approximation  theory  for  infinite 
dimensional  discrete-time  LQG  problems  to  the  problem  here  is  discussed, 
and  numerical  examples  are  presented. 
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Computational  Methods  for  the  Solution  of  Infinite  Dimensional  Discrete-Time 
Regulator  Problems  with  Unbounded  Input,  I.  G.  Rosen  and  M.  A.  Lie. 

An  approximation  framework  for  the  closed-loop  solution  of 
discrete-time  linear-quadratic  regulator  problems  for  infinite  dimensional 
systems  with  unbounded  control  inputs  is  developed.  Sufficient  conditions 
for  the  convergence  of  approximations  to  Riccati  operators  and  feedback 
gains  which  characterize  the  optimal  control  law  are  provided.  General 
theories  for  abstract  partial  differential  systems  with  boundary  control  and 
distributed  systems  with  control  delays  are  developed.  Spline-based 
schemes  and  numerical  results  for  heat  and  beam  equations  with  boundary 
control  and  a  hereditary  system  with  delayed  control  are  presented  and 
discussed. 

Approximation  of  Discrete-Time  LQR  Problems  for  Boundary  Control  Systems  with 
Control  Delays,  I.  G.  Rosen. 

In  this  short  note  we  consider  the  extension  and  application  of  the 
approximation  theory  for  discrete-time  linear-quadratic  regulator  problems 
with  either  bounded  or  unbounded  inputs  we  developed  earlier  to  boundary 
control  systems  with  control  delays.  We  synthesize  our  earlier,  existing 
results  for  distributed  systems  with  boundary  controls  and  for  systems  with 
control  delays  into  a  theory  which  is  applicable  to  systems  that 
simultaneously  exhibit  both  forms  of  unbound^  input. 
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ABSTRACT 


We  consider  the  approximation  of  optimal  discrete-time  linear  quadratic  Gaussian  (LQG) 
compensators  for  distributed  parameter  control  systems  with  boundary  input  and  unbounded 
measurement.  Our  approach  applies  to  problems  tliat  can  be  fonnulated  in  a  state  space  on  which 
both  the  discrete-time  input  and  output  operators  are  continuous.  Approximating  compensators  are 
obtained  via  application  of  the  LQG  theory  and  associated  approximation  results  for  infinite 
dimensional  discrete-time  control  systems  with  bounded  input  and  output.  Numerical  results  for 
spline  and  modal  approximation  schemes  used  to  compute  optimal  compensators  for  a  one 
dimensional  heat  equation  with  either  Neumann  or  Dirichlet  boundary  control  and  pointwise 
measurement  of  temperature  are  presented  and  discussed. 
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In  this  paper  we  develop  an  approximation  theory  for  the  computation  of  optimal  discrete-time 
linear  quadratic  Gaussian  (LQG)  compensators  (combined  feedback  control  law  and  state  estimator) 
for  distributed  parameter  systems  with  boundary  input  or  control  and  unbounded  output  or 
measurement.  In  a  continuous  time  setting,  boundary  input  typically  results  in  an  unbounded  input 
operator.  That  is,  the  system's  input  operator  maps  the  control  input  into  a  space  larger  than  the 
state  space  in  which  the  open-loop  system  is  usually  formulated.  In  the  discrete-time  case,  on  the 
other  hand,  for  a  wide  class  of  distributed  systems,  the  resulting  input  operator  is  bounded  on  the 
usual  underlying  state  space.  An  unbounded  output,  or  measurement,  operator  has  domain  smaller 
than  the  usual  open-loop  state  space. 

For  continuous  time  systems,  Pritchard  and  Salamon  (1987)  have  established  an  abstract 
semigroup  theoretic  framework  for  treating  the  linear  quadratic  regulator  problem  (control  only)  for 
infinite  dimensional  systems  with  unbounded  input  and  output  operators.  Their  approach  is  based 
upon  a  weak  or  distributional  formulation  of  the  Riccati  equations  which  characterize  the  optimal 
feedback  control  laws  in  an  appropriate  dual  space .  Curtain  (1984)  provides  a  procedure  for  the 
design  of  finite  dimensional  compensators  for  parabolic  systems  with  unbounded  control  and 
observation.  In  (Curtain  and  Salamon,  1986)  a  finite  dimensional  compensator  design  procedure 
for  a  wider  class  of  infinite  dimensional  systems  with  unbounded  input  (but  bounded  output) 
including  hereditary  systems  with  control  delays  and  partial  differential  systems  with  boundary 
control  is  developed.  Lasiecka  and  Triggiani  have  looked  at  linear  regulator  problems  for  parabolic 
(1983a,  1987a)  and  hyperbolic  (1983b,  1986)  systems  with  boundary  control  and  obtained, 
among  other  things,  global  and  local  regularity  results  for  the  optimal  controls  and  state 
trajectories.  In  (Lasiecka  and  Triggiani,  1987b)  Galerkin  approximations  and  an  associated 
convergence  theory  for  closed-loop  solutions  to  regulator  problems  for  parabolic  systems  with 
Dirichlet  boundary  input  are  studied.  A  more  complete  survey  of  the  boundary  control  literature 
including  references  to  some  of  the  poineering  work  in  this  area  can  be  found  in  (Pritchard  and 
Salamon,  1987). 


1 


In  our  treatment  here,  we  consider  the  discrete-time  problem  (i.e.  piecewise  constant  input  and 
sampled  output).  Our  interest  in  the  discrete-time  or  digital  fomiulation  is  motivated  by  1)  the  fact 


that  it  represents  a  more  accurate  or  realistic  description  of  how  the  linear-quadratic  theory  for 
distributed  systems  would  actually  be  applied  in  practice,  and  by  2)  how  the  boundedness  of  the 
discrete-time  input  operator  in  the  usual  underlying  state  space  facilitates  the  development  of  an 
approximation  theory  which  can  simultaneously  handle  both  unbounded  input  and  unbounded 
output.  Our  approach  is  based  upon  an  application  of  the  theory  we  developed  earlier  in  (Gibson 
and  Rosen,  1985  and  1986)  for  the  approximation  of  optimal  discrete-time  LQG  compensators  for 
infinite  dimensional  systems  with  bounded  input  and  output  Our  results  are  applicable  to 
boundary  control  systems  in  which  a  restriction  of  the  state  transition  operator  and  the  discrete-time 
input  operator  are  bounded  on  a  space  on  which  the  output  operator  is  bounded  as  well.  To 
illustrate  our  approach,  in  this  paper  we  describe  in  detail  the  application  of  our  theory  to  a  one 
dimensional  heat  equation  with  either  Neumann  or  Dirichlet  boundary  control  and  pointwise 
measurement  of  temperature.  Elsewhere  (see  Gibson  and  Rosen,  1987)  we  have  applied  our 
results  to  develop  approximation  schemes  for  the  computation  of  optimal  LQG  compensators  for 
flexible  structures  (i.e.  Euler-Bemoulli  beams)  with  shear  force  input  at  the  boundary  and  a 
pointwise  measurement  of  strain. 

An  outline  of  the  remainder  of  the  paper  is  as  follows.  In  Section  2  we  describe  an  abstract 
framework  for  the  study  of  boundary  control  systems  and  their  discrete-time  formulation.  In 
Section  3  we  review  the  LQG  theory  for  infinite  dimensional  discrete-time  systems  and  associated 
abstract  approximation  results.  In  the  fourth  section,  we  discuss  spline  and  modal  subspace  based 
approximation  schemes  for  the  heat  equation  example.  Section  5  contains  some  concluding 
remarks. 

2.  The  Boundary  Control  System  and  its  Discrete-Time  Formulation 

We  employ  a  semigroup  theoretic  formulation  that  has  been  used  previously  for  a  class  of 
abstract  boundary  control  systems.  See,  for  example,  (Curtain  and  Salamon,  1986).  Let  W,V  and 
H  be  Hilbert  spaces  with  W  and  V  densely  and  continuously  embedded  in  H.  We 


consider  boundary  control  systems  of  the  fomi 


(2.1) 

w(t)  =  Aw(t), 

t>0 

(2.2) 

o 

II 

o 

(2.3) 

Fw(t)  =  v(t). 

t>0 

(2.4) 

y(t)  =  Cw(0. 

t>0 

where  A  e  33  (W.H),  the  boundary  input  operator  F  is  an  element  in  33(W,R"')  and  the  output 
operator  C  is  an  element  in  33(V,RP).  Note  that  the  operator  A  need  not  be  the  Laplacian.  Our 
choice  of  A  to  denote  a  general,  most  often  differential,  operator  satisfying  the  conditions  set  forth 
below  is  consistent  with  the  notation  used  in  earlier  treatments  of  boundary  control  systems 
elsewhere  in  the  literature. 

We  assume  that  1)  F  is  suijective  and  its  null  space,  Tl(F)  =  {cp  e  W:  Ftp  =  0},  is  dense  in 
H,  2)  the  operator  Cl,  defined  to  be  the  restriction  of  the  operator  A  to  TI>(F),  is  a  closed  operator 
on  H  and  has  non-empty  resolvent  set  and  3)  for  each  T  >  0,  all  Wq  e  W,  and  v  £  C’(0,T;  R^) 
with  Fwq  =  v(0),  there  exists  a  unique  w  e  C([0,T];  W)  H  CH[0,T]:  H)  which  depends 
continuously  on  Wq  and  v  and  which  satisfies  (2.1)  -  (2.3)  for  each  t  e  [O.T].  It  then  follows  (see 
Hille  and  Phillips,  1957)  that  the  operator  C( :  Dom  (Cl)  c  H  — >  H  given  by  Cltp  =  Acp  for  cp  E 
Dom(Cl)  =  71(0  is  the  infinitesimal  generator  of  a  Gq  semigroup,  {ir(t):  t>  0},  of  bounded 
linear  operators  on  H. 

Define  the  space  Z  as  the  dual  of  Dom(Cl  )  where  the  norm  on  Dom  (Cl  )  is  taken  to  be  the 
graph  Hilbert  space  norm  associated  with  the  operator  Cl*.  Then  H  is  densely  and  continuously 
embedded  in  Z  and  { tr(t) :  t  ^  0}  can  be  uniquely  extended  to  a  Gq  semigroup  of  bounded  linear 
operators  on  Z.  Its  generator  is  the  extension  of  the  operator  Ci  to  an  the  operator  ^  in  G(H,Z) 
given  by  (  Q  (p)(v)  =  <(p,  Cl  for  <p  €  H  and  e  Dom(Ci  ). 

Since  F  was  assumed  to  be  a  suijcction,  it  has  a  right  inverse.  Let  F^  :  R'”  — >  W  be  any  right 
inverse  of  F .  Since  Dom  ( F^)  =  R'",  we  have  F^  £  G(R'^,  W).  For  v  £  R'",  we  define 
73  £  33 (R*^,  Z)  by  Bv  =  (A  -  ^)  F^  v.  If  F^  and  F^  are  two  distinct  right  inverses  of  F  then 
B.(  F+  -  F+ )  c  71(F).  Since  ^  coincides  with  A  on  71.  (F) ,  it  follows  that  the  operator  73  is 


well  defined.  It  can  be  shown  (see  Curtain  and  Salamon,  1986)  that  for  each  Wq  e  fi  and  v  e 
1^(0, T;  R'")  there  exists  a  unique  w  e  C([0,T);  H)  H  (0,T;  Z)  which  depends  continuously  on 
Wq  and  V  and  which  satisfies 

w(t)  =  Ctw(t)  +  (I3v(t),  t  >  0 
w(0)  =  Wq 

in  Z.  The  function  w  is  given  by 


(2.5)  w(t)  =  tT (t)Wg  +  J  cr(t  -  s)  Bv(s)ds,  t  >  0 


and  is  referred  to  as  a  weak  solution  to  the  boundary  control  system  (2.1)  -  (2.3). 

The  discrete-time  formulation  of  (2.1)  -  (2.4)  is  found  by  considering  piece  wise-constant 
controls  of  the  form 


(2.6)  v(t)  =  Ujj,  t  e  [kt,  (k -t- l)x),  k  =  0,1,2,... 

where  t  denotes  the  length  of  the  sampling  interval.  Let  Wj^  =  w(kx) ,  k  =  0,1,2,...  where  w(  ) 
is  the  unique  weak  solution  to  (2.1)  -  (2.3)  given  by  (2.5)  corresponding  to  Wq  e  H  and  input  v 
given  by  (2.6).  (We  note  that  with  piecewise  constant  input  of  the  form  (2.6),  the  solution  w  is  in 
fact  a  strong  solution  on  each  subinterval  [kx,  (k+l)x].)  For  each  k  =  0,1,2,...  define 
e  C([kx,  (k  +  l)x];  H)  by  z^(t)  =  w(t)  -  F^u^,  t  e  [kx,  (k  +  l)x].  Then 

zic(0  =  '''(t)  =  Ciw(t)  +  ^3u^ 

=  Cl  2,^(0  -f-  (Cl  + 

=  Ciz^(t)  +  AP-u,^,  te(kx,  (k+ l)x], 
z^(kx)=  w^-p-u^. 

Therefore 


V-  ••  V, v.v.‘ 


'Vi  =^((^+i)^)  +  r^“k 

=  3r(t)(w^  -  p-uj.)  +  r  3^(s)  A  r+u,^ds  +  r+u^ 


=  3r(x)w^  +  (I-0'(t))r+u^+  1^  a'(s)APu^ds. 


'^k+i  +  k  =  0,1,2,... 


WqE  H 


where  T  e  a3(H)  and  B  e  aSCR"*,  H)  are  given  by  T  =  ^(x)  and  B  =  (I  -  iT (x))  1^  + 


tT  (s)An*'ds  respectively. 


We  note  that  as  in  the  case  of  the  continuous-time  input  operator  the  discrete-time  input 


operator  B  is  well  defined  and  does  not  depend  upon  a  particular  choice  for  P.  Indeed  if  Bj  and 
B2  are  the  input  operators  which  correspond  to  the  choices  and  Fj"*"  then  for  u  e  R*"  we  have 


(Bj  -  B2)u  =  (I  -  T(x))(  F^  -  r2+)u  -t-  r  0'(s)A  (Fi^  -  F2+)uds. 


But  (F'*’  -  F^  )u  e  U(F)  =  Dom(C()  and  therefore 
1  2 


T  T 

j  or  (s)A(P  -  P  )uds  =  J  U'(s)  (U(P  -  P  )  uds 


=  ^(s)(r;  -  q  )uds  =  (O'(x)  -  ixq  - q )  u . 


In  addition,  if  P  is  chosen  so  that  in(P)  c  %(A),  B  takes  on  the  particularly  simple  form  B  = 


(I  -  0^ (x))P .  It  is  worth  noting  that  a  simple  calculation  reveals  that 


B  =  /  ^  0^(5)  Bds 


in  agreement  with  the  standard  technique  for  obtaining  the  discrete  or  sampled  time  formulation  of 


a  continuous  time  system  in  either  a  finite  dimensional  or  bounded  input  setting. 


It  is  our  intention  here  to  apply  the  approximation  theory  we  developed  earlier  in  (Gibson  and 
Rosen,  1986)  for  the  design  of  optimal  discrete-time  LQG  compensators  for  infinite  dimensional 
systems  with  bounded  input  and  output  operators.  We  therefore  require  the  additional  assumptions 
that  4)  T  =  cr(x)  e  33(V)  and  5)  (fi.(P')  c  V.  Although  not  all  boundary  control  systems  we 
might  formulate  would  satisfy  these  conditions,  there  are  many  interesting  and  important  systems 
which  do  (see,  for  example.  Section  4  below  and  Gibson  and  Rosen,  1987).  In  this  case,  the 
control  system  (2.1)  -  (2.4)  takes  the  form 

(2.7)  =  Tw^  Bu^.,  k  =  0, 1 ,2,... 

(2.8)  WoeV 

(2.9)  yk  =  Cwi.,  k  =0,1,2,...  . 


3.  LOG  Theory  for  Infinite  Dimensional  Discrete-Time  Systems  and  Finite 

Dimensional  Approximation 


3.1.  The  Infinite  Dimensional  Probelm 

The  discrete-time  linear-quadratic  regulator  problem  for  the  boundary  control  system  (2.1)  - 
(2.3)  is: 


Findu  =(u\)k^oe  ooj  R"’)  which  minimizes  the  quadratic  performance  index 


J(u)= 

fco 


where  Qe  33  (V)  is  self-adjoint  and  nonnegative,  R  is  a  symmetric  positive  definite  mxm  matrix 


and  the  state  w={wt}~  evolves  according  to  the  recurrence  (2.7),  (2.8). 
“  x-o 


An  optimal  control  exists  for  each  initial  condition  Wq  if  and  only  if  the  operator  algebraic 
Riccati  equation 


(3.1)  n=  T*(n-nB(R  +  B*nB)'B*n)T  +  Q 


has  a  bounded  nonnegative,  self-adjoint  solution  n.  In  this  case,  the  optimal  control  has  the 

feedback  form  Uj.  = -Fwjj  where  F  =  (R  +  B  nB)'*B  FIT.  A  control  (sequence')  u  is  admissible 

for  the  initial  condition  Wq  if  the  corresponding  J(u)  is  finite.  If  there  exists  an  admissible  control 

for  each  initial  condition,  then  (3.1)  has  a  bounded  nonnegative,  self-adjoint  solution.  If  each 

admissible  control  for  each  initial  condition  drives  the  state  to  zero  asymptotically,  then  there  exists 

at  most  one  bounded  nonnegative,  self-adjoint  solution  to  (3.1).  The  optimal  trajectory  w''  = 

{w*,, }  evolves  according  to  w*  =  k  =  0,1,2,...,  where  the  closed  loop  state  transition 
k=0  k 

operator  S  e  JS(V)  is  S  =  T  -  BF.  If  Q  is  coercive,  then  S  has  spectral  radius  less  than  one  and  is 
uniformly  exponentially  stable.  From  the  finite  dimensionality  of  the  control  space  we  obtain 


Ut  =  -<f,w!>v,  k  =  0,1,2,... 


where  f  =  (fpf2,...,fn^)*^  e  x  V  is  called  the  optimal  functional  feedback  control  gain. 

j=l 

The  results  stated  here  for  the  optimal  linear-quadratic  regulator  problem  are  summarized  from 
(Gibson  and  Rosen,  1985). 

When  only  a  finite  dimensional  measurement  y  =  { ^  =  o  of  infinite  dimensional  state  w 
is  available  (recall  (2.9)),  a  state  estimator  or  observer  is  required.  For  a  given  input  sequence  u 
and  corresponding  output  sequence  y,  the  optimal  LQG  estimator  is 


(3.3)  w^+l  =  Tw^-^Bu^-t-F(y^^-Cw^^},  k  =  0,1,2,... 

(3.4)  Wq  e  V 

where  the  optimal  estimator  or  observer  gain  F  e  JC(RP,V)  isF  =  TITC  (R-^C^C  )"' with 
n  6  3^(V)  the  minimal,  self-adjoint,  nonnegative  solution  (if  one  exists)  to  the  operator  algebraic 
Riccati  equation 


^  /N  .  ^  ^ 


n  =  T(  n  -  n  c*(r  -t-  c  n  c*)  ’c  n  )t*  -t-  q 


Since  F  6  (RP,V),  it  has  the  representation 


Fy  =  f^y  ,  y  E  RP 


where  f  =  (f  j,  f  2.--.  fp)  ^  e  x  V  is  called  the  optimal  functional  observer  gain. 

j=i 

The  operator  Q  e  I5(V)  is  self-adjoint,  nonnegative  and  the  pxp  matrix  R  is  symmetric, 
positive  definite. 

In  a  stochastic  setting,  the  operator  Q  and  the  matrix  R  are,  respectively,  the  covariance 
operator  and  covariance  matrix  for  uncorrelated,  zero-mean,  stationaiy,  Gaussian  white  noise 
processes  that  force  the  state  and  corrupt  the  measurement.  In  this  case,  if  Q  is  trace  class,  (3.3), 
(3.4)  is  the  infinite  dimensional  analog  of  the  discrete-time  Kalman-Bucy  filter.  In  a  strictly 
deterministic  setting,  Q  and  R  are  assumed  to  be  determined  via  engineering  design  criteria 
such  as  stability  margins,  robustness  of  the  closed-loop  system,  etc. 

Replacing  operators  in  the  control  problem  with  the  adjoints  of  the  appropriate  operators  in  the 
estimator  problem  yields  the  usual  duality  between  the  LQG  optimal  control  and  estimator 
problems.  Hence  sufficient  conditions  for  existence  and  uniqueness  of  solutions  of  (3.5)  and  the 
closed-loop  estimator  stability  properties  are  analogous  to  the  results  for  the  control  problem.  In 
particular,  if  Cj^  =  wj^  -  Wj^,  then  Cq,  k  =  0, 1,2, ...,  where  S  =  T  -  FC,  and  a  sufficient 

condition  for  5  to  be  uniformly  exponentially  stable  is  that  0  be  coercive. 

The  optimal  LQG  compensator  consists  of  the  state  estimator  in  (3.3)  and  (3.4)  and  the  control 


Ut*  =  -F  w  *  ,  k  =  0,1,2,. 


The  resulting  closed-loop  system  is  given  by 


*11;.=  k  =0,1,2,... 


where  *11/^  =  (w^.,  w*^)^  with  { ,, o  trajectory  that  results  from  the  input  (3.6) 

and  ^  e  2;(VxV)  is 


\rm  wn  v-  . 


FC  T-BF-FC 


It  is  easy  to  show  that  tlie  spectrum  of  ^  is  given  by  a(/6)  =  a(S)  U  a(^),  so  that  the  stability  of 
the  closed-loop  plant-compensator  system  is  determined  by  the  stability  of  the  plant  with  full  state 
feedback  and  the  stability  of  the  estimator  error. 


For  each  N  =  1,2,...,  let  Vjyf  be  a  finite  dimensional  subspace  of  V  and  let  be  a  bounded 
linear  mapping  from  V  onto  Vj^  (for  example,  the  orthogonal  projection  with  respect  to  either  the  V 
or  H  inner  product).  Let  Tj^t,  Qj^,  Qi,j  e  S(Vj^),  Bj^e  3S(R'”,Vj^)  and  Cj^  e  JS(Vj^,RP)  and  set 

F„=  (R  +  B*n„B^,)Vn„TN 

and 

^N'T^nNCjCR-i-c^fiNCV' 

where  ITj^  and  ITj^  are  the  minimal,  self-adjoint,  nonnegative  solutions  (assuming  that  they  exist)  to 
the  finite  dimensional  operator  algebraic  Riccad  equations 

(3.7)  rij^  =  ITj^ '  ITf^Bf^R  -H  B^rifjBj^)  *  B^FIj^jTj^  +  Qj^ 
and 

(3.8)  On  =  T^C  -  fl^  C‘  (R  +  C^flN  C"  )-'Cn  )T*+ 

respectively.  The  approximating  optimal  compensator  is  given  by 


G*  =-Fnw"  .k  =  0,1,2,.. 
n;c  n,x 


9 


where  w*  =  { w  }~  is  detemiined  according  to  the  approximating  observer 
n  n;ix  =  o 


w*  =  Tf.  w*  +Bj^u*  +  Fw{y*  -  Cw  w*  },  k;  =  0,1,2. 
NX  +  i  N  n;c  N^x  N,x  ^  n;c 

'^N,0  ^  '^0  £  • 

The  measurements  y  are  given  by  y  =  Cw^.^  ,  k  =  0,1,2,...  where 
N  N ,  X  ' 


WN.kH=Twj^.k-Bu*  ^  .  k  =  0,1,2,... 


Wn.o  =  Wq  . 


The  resulting  closed-loop  system  is  given  by  “^nq*  ^  =  0,1,2,...  where 


‘'^N.k  =  ('^N.k-  V  and  e  (VxV^)  is  given  by 

*  *  N  Jt 


/S  ^ 

Let  Sj,^  =  Tj^  -  and  Sj,^  =  Tj^  -  Fj^Cj^  and  assume  that  Pj,^  ->  I  strongly  on  V  as 
N  ->  oo.  Assume  further  that  Tj^Pj^  —>  T,  T*  Pj^  — >  T*,  QjsjPj^  — »  Q  a!nd  Qj,|P[^  ->  Q 
strongly  on  V  and  that  Bj,^  B  and  Cj^Pj^  ->  C  in  norm  as  N  ->  «».  If  the  pairs  (Tj^,  Bj^) 
and  (T*  ,  C* )  are  uniformly  exponentially  stabilizable  and  the  pairs  (Tj^,  Qj^)  and  (T*  .  Q^)  are 
detectable  (see  Kwakemaak  and  Sivan,  1972)  then  there  exist  unique,  self-adjoint,  nonnegative 
solutions  rij^^and  to  the  algebraic  Riccati  equations  (3.1)  and  (3.5).  If  ITj^  and  flj^  are 
bounded  from  above  uniformly  in  N,  then  flj^Pj^i  and  Hj^^Pn  converge  weakly  to  IT  and  If, 
respectively,  as  N  — >  <»  , 

If,  in  addition,  Sj,j  and  Sj^  are  uniformly  exponentially  stable,  uniformly  with  respect  to  N, 
then  rij^P}^  and  converge  strongly.  Weak  convergence  of  Hj^^Pj^  to  FI  yields  strong 

convergence  of  Fj,jPj,j  to  F  and  Sj^Pj^  to  S.  If  flj^^Pj^  converges  strongly  then  Fj^Pj.^  F  in 


iC*. 

i 


'‘.V, 


nonn.  Weak  convergence  of  to  IT  yields  weak  convergence  of  Fj^  to  F  and  to 

S.  When  FIj^jPj^  — >  fl  strongly,  then  Fj^  — >  F  in  norm  and  Sj^Pj,j  — >  S  strongly  in  V  as 
N ->  oo.  Finally,  if  is  the  mapping  of  VxV  onto  VxVj^  given  by  ?°j^(Wj.  W2)  =  (Wj, 

PnW2),  then  ^J^PI^J  ->  FI  weakly  or  strongly  is  sufficient  to  conclude  that  A  weakly 

or  strongly  depending  only  upon  whether  flj^Pj^  ->  fl  weakly  or  strongly  as  N  ^  .  Under 

appropriate  additional  hypotheses  on  the  spectral  properties  of  the  open-loop  system  and  on  the 
approximation  scheme,  it  is  possible  to  show  that  converges  to  ,6  in  norm.  (We  have  been 

able  to  obtain  such  a  result  only  for  modal  approximations.)  Norm  convergence  of  the  closed-loop 
state  transition  operators  is  sufficient  to  conclude  that  uniform  exponential  stability  of  /6  implies 
uniform  exponential  stability  of  for  all  N  sufficiently  large  (see  Gibson  and  Rosen  1986). 

In  practice,  the  finite  dimensional  approximating  subspaces  are  often  constructed  using  any 
of  a  number  of  common  finite  element  bases,  e.g.  polynomial  and  hermite  spline  functions,  mode 
shapes,  orthogonal  polynomials,  etc.  For  the  discrete-time  boundary  control  systems  of  interest 
to  us  here,  the  approximations  to  T  and  B,  and  are  obtained  by  approximating  the 
continuous  time  semigroup,  { 3^(t) :  t  ^  0),  by  a  semigroup  of  bounded  linear  operators  on 
{T j^(t):  t  ^  0).  In  fact  it  is  the  infinitesimal  generator  (!l  of  the  semigroup  {O' (t):  t  >  0}that  is 
approximated  by  a  bounded  linear  operator  on  Vj^  with  ( «rj^(t):  t  >  0}  then  being  defined  by 

Tj^ft)  =  exp  (Cl^t).  t  >  0.  With  =  0'j^(x)  and  =  (I  -  0'j,^(T))Pj,jr+  -t-  0'j^(s)PNAr^ds,  the 
required  convergence  can  usually  be  proved  using  the  Trotter-Kato  semigroup  approximation 
result  (see  (Kato,  1966)  and  (Pazy,  1983) ).  The  approximations  to  Q,  Q  and  C,  and 

Cj^ ,  respectively,  typically  are  taken  to  be  Qj^  =  P^^Q,  Qj^  =  Pj^Q  and  =  CPj^ . 

denote  a  basis  for  Vj^  and  set  0^  =  ((pN  ,(p^,...,  9^  )^e  \  Vj^. 

)  J  ™  I  1  2  N  j_j 

Adopting  the  convention  that  [L]  denotes  the  matrix  representation  with  respect  to  the  basis 
(9!^}  j*'  lor  a  linear  operator  L  with  domain  and/or  range  in  Vj.^,  we  find  that 
[F^,]  =  (R  -H  [B^f  0  N  [B^D-HB^f  0  and  [F^,]  =  [T^]  (R  [C^] 

0  '  where  0^and  0^  are  the  unique,  synunetric,  nonnegative  solutions  to  the 

’’n  ^  rnatrix  algebraic  Riccati  equations 


(3.10) 


0^  =  [TnF  (0N  -  0N[B^](R  +  0  N[BN])-l[BNf  0N)ITn] 

+  MN[Q^] 


(3.11) 


0N  =  [TnK  0^-  0^[CNf  (R  +  [CN]0^[CNf)-‘[CN]0^)[TN]T 
+  [Qn](M^)-‘  . 


The  matrix  is  the  Pj^  x  Pj^  Gram  matrix  ,  (0^)^>Y  , 


If  =  W*_^e  R^N.then  >  -[F^]  ,  k  =  0,1,2,...  with 

.r  "  '®n)  [FnKJ;,  -  [C«)  w;_^),  k  =  0,1,2,... 


w;;^^=(M'')-<«i>i',w„>v . 


Ill 

The  approximatipg  optimal  fupctiopal  feedback  control  gain,  =  (f^,  e  x 

j=i 

are  given  by  f^=  [Fj^](M^)'^<I>^  and  the  approximating  optimal  functional  observer  gain 


/s  ^  A 


f^=(f^,  f^,...,  f^^  e  X  Vj,^  by  f  =  [Fj^l'^O^,  Ifnj,^Pj,^  ->n  weakly  (strongly) 

j=i 

At  A 

then  f^  ->  fj,  i=  1,2,...,  m  weakly  (strongly)  in  V.  If  -^11  weakly  (strongly)  then 

f^  — >  fj  ,i  =  1,2,...,  p  weakly  (strongly)  in  V.  If  the  injection  V  c  H  is  compact,  then  fj^  fj , 

/V  ^  A 

i=l,2 . mand  fj  -»  fj.i  =  l,2,...,p  strongly  in  H  if  rTf^Pj,!  and  Hj^P^  converge  only 

weakly. 


We  consider  the  one-dimensional  heat  equation 
3\v  3  w 

(4.1)  - (t,x)=  a — -(t,x),  0<x<l, 


I 


where  a  >  0,  with  the  homogeneous  Dirichlet  boundary  condition 


(4.2)  w(t,0)  =0,  t  >  0. 


and  either  the  Neumann  boundary  control 


(4.3)  _(t.l)  =  v(t).  t>0. 

dx. 

or  the  Dirichlet  boundary  control 


(4.4)  w(t.l)=v(t).  t>0. 

where  v  e  ^(O, «).  For  output  we  take  a  temperature  measurement 

(4.5)  y(t)  =  w(t.O,  t>0. 

at  some  fixed  point  C  ^  (0,  1).  Initial  conditions  for  these  systems  have  the  form 

(4.6)  w(0,  x)  =  Wq  (x),  0  ^  X  ^  1 

where  Wq  e  L2  (0, 1). 

Although  the  two  control  systems  above  appear  to  be  similar,  they  are,  in  fact,  quite  different  and 
must  be  treated  separately.  We  begin  with  the  more  straight  forward  of  the  two,  Neumann  boundary 
control.  Let  H  =  Lj  (0, 1),  V  =  HJ^(0,  1)  =  {<p  E  H^O,  1) :  9  (0)  =  0}  and 
W  =  h2  (0, 1)  n  (0,  I).  With  H  endowed  with  the  usual  inner  product,  V  with  the  inner 

product  <  9,  V  >v  ~  -^0  W  with  the  inner  product  <  9,  \|/  ^  Jg  Di  DJ  9,  we 

j=l 

have  the  continuous  and  dense  embeddings  WcVcHcV’c  W.  Define  A  e  2(W,  H),  F  e  33  (W, 
R’)  by  a  9  =  a  D^9,  r9  =  D9(l)  and  C9  =  9(Q  respectively.  With  these  definitions  the  boundary 


IHJPJiLIPLILIlim 


control  system  (4.1)  -  (4.3),  (4.5),  (4.6)  has  the  fomi  (2.1)  -  (2.4).  The  operator  Ct:  Dom  (Q)  c 

H— >H  is  given  by  Ct.(p  =  a  D^tp  for  tp  e  {(p  e  H^(0,  1);  (p(0)  =  D(p(l))  =  O).  It  is  densely  defined, 

negative  definite,  self-adjoint  and  it  is  the  infintesimal  generator  of  a  unifomily  e.xponenlially  stable 

analytic  semigroup  {(T  (t):  t  >  0}  of  bounded,  self-adjoint  linear  operators  on  H.  Also,  {O' (t):  t  > 

0)  is  a  unifomily  exponentially  stable,  analytic  semigroup  of  bounded,  self-adjoint  operators  on  V 

with  generator  ^  given  by  ^cp  =  Cttp  for  tp  e  {(p  e  (0,  1):  (p(0)  =  D(p(l)  =  D^(p(0)  =  0}. 

Choosing  e  2S(R^  W)  as  (  P''u)(x)  =  xu  for  x  e  [0,  1],  we  have  R(  P*”)  c  V,  7\.(  P*")  c  Tl(A) 

and  that  conditions  1)  -5)  given  in  Section  2  are  satisfied.  For  the  optimal  control  and  estimator 

problems,  we  take  Q  =  ql,  0  =  I.  R  =  r  and  ft.  =  f  where  I  is  the  identity  on  V,  q,  q  >0  and  r,  f 

>  0.  The  uniform  exponential  stability  of  the  semigroup  {ir(t):  t  >  0}  on  V  implies  that  the  algebrai 

Riccati  equations  (3.1)  and  (3.5)  admit  unique  bounded,  nonnegative,  self-adjoint  solutions  n  and  fl 

respectively.  The  optimal  control  (3.2)  takes  the  form 

1 

(4.7)  u*  =  -  j  DfDw*.  k=  0,1,2,... 

0 


where  the  optimal  functional  feedback  control  gain  f  and  the  optimal  functional  observer  gain 

A 

f  are  elements  in  (0,1). 

We  construct  an  approximation  scheme  using  a  linear  spline  based  Ritz-Galerkin  approach.  For 

each  N  =  1,2,... ,  ^  denotes  the  usual  linear  spline  or  "hat"  functions  defined  on  the  interval 

[0,1]  with  respect  to  the  uniform  mesh  (0,  1/N,  2/N . 1 }.  We  discard  the  element  centered 

at  X  =  0,  tp^,  set  Vf^  =  span  {(p|^}  ^  and  choose  to  be  the  orthogonal  projection  of  V  onto  Vj^ 

with  respect  to  the  V  inner  product  Hence  Vj,,  is  an  N  dimensional  subspace  of  V. 

For  (p  e  Dom  (Ct),  >  altply  >  319!^^  and  therefore  0  E  p((U)  and  H  — >  Dom(Ct) 

satisfies  I  Cl' ^9Iy  ^  We  define  Clj.,:Vf^  — >  Vj^  as  the  inverse  of  the  operator 

Cp^  =  PjsjCp'  restricted  to  Vj^.  The  operator -Cl"  *  is  positive  definite  because 
N  N 

<  Cl*  ^  9j.^,  9n>v  =  for  9^  e  V^,  and  it  is  self-adjoint  since  <  Cl"  ‘  9^^,  Vn>v  = 

M  H  H 

<Pf^Cl'' 9j,j,  =  <Cl'*  9fj,  V|rj^>V  =  a'*  <  9n’ •  Hence  the  operator  Cl^  is  well  defined 

and  self-adjoint  For  9^^  e  and  =  Clj^9j^  ,  the  estimate 


^  -aia-VN'v^  -^'Pn^'Wn'y  = 

=  -a|cpN'v 


implies  that  Clj^  is  the  infinitesimal  generator  of  a  Cq  semigroup  { iTj^Ct)  :  t  >  0)  of  bounded, 
self-adjoint  linear  operators  on  Vj,^  satisfying  I  cTj^ft)!  <  e'^  ^  ,  t  >  0. 

It  can  be  shown  that  a<(p  ,  \ir>y  =  <(-Ct)*^^q),  (-Ct)^^Y>pj .  It  then  follows  that  the  matrix 
representation  for  the  operator  Ctj^  with  respect  to  the  basis  {tp!^}  is  [Ctf,j  ]  = 

-a<  (p^,  (p!^>H  <  tpf^,  tpf^>v  •  This  agrees  with  the  system  matrix  derived  by  a  standard 
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Ritz-Galerkin  finite  element  approach.  Note  that  even  though  C|j^  is  defined  to  be  the  inverse  of  the 
operator  restricted  to  the  space  Vjv^,  computing  its  matrix  representation  does  not  require 

either  or  Ci'^  explicitly.  In  general,  the  same  approach  can  be  used  to  obtain  an  operator 
representation  for  the  Ritz  Galerkin  approximation  to  any  self-adjoint  coercive  operator. 

Let  denote  the  interpolation  operator  from  V  onto  defined  by  (dj^  (p)(j/N)  =  (p(j/N),  j  =  1 , 
2,...,  N .  Then  for  9  e  W,  elementary  approximation  properties  of  linear  interpolatory  spline 
functions  (see  (Schultz,  1971))  imply 


I(Pn-D9'v  ^ 

Nk 

and  therefore,  since  W  is  dense  in  V,  that  strongly  on  V  as  N— >  .  Also,  it  follows  that 

Cl^=  Pj^Cl^— >  strongly  on  V  as  N—>  «>  .  If  we  define  Tj^  =  O' j^(t),  then  the  Trotter-Kato 
approximation  theorem  yields  that  Tj^Pj^  ->  T  strongly  on  V  as  N-> and,  since  T*  =  T  =  Tft) 
and  TJI^  =  Tj^  =  O' j^t)  that  Tj!^Pj^  ->  T*  strongly  on  V  as  N->  <». 

Since  (r\,(r^  )  c  Vj,j  (recall  that  (r^u)(x)  =  ux,  0  ^  x  ^  1),  we  define  the  approximating  input 
operators  by  =  (I  -  O' ^(x))!^  and  set  Q,^  =  ql,  =  q  I  and  Cj^  =  C.  The  strong  convergence 
of  Pf^  to  the  identity  and  Tj^^Pj^  to  T  together  with  the  finite  dimensionality  of  the  domain  of  B  and  the 
range  of  C  are  sufficient  to  conclude  that  Q^P^  — >  Q,  QnPjsj  Q  strongly  on  V  and  that  Bj,^  B 
and  — >  C  in  norm  as  N  ->  <». 

The  uniform  exponential  stability  of  the  semigroups  {  O',.!!)  ;  t  >  0)  implies 


(4.8)  IT|^iY  =  l(r  )\<r>‘.  k  =  0.1,2.... 


with  r  =  e  <  1.  Consequently  the  pairs  (Tj.^.Bj^)  and  (Tj.^,  C*  )  are  unifonnly  exponentially 

A 

stabilizable  and  the  pairs  (Tj.j,Qj.j)  and  (Tj^,  Qj.^)  are  detectable.  It  follows  that  there  exist  unique 

-iN 

self-adjoint,  nonnegative  solutions  FIi^  and  n|.j  to  the  finite  dimensional  algebraic  Riccati 
equations  (3.7)  and  (3.8)  respectively.  The  uniform  exponential  bound  (4.8)  with  r  <  1  imples  that 

A>. 

the  zero  control  yields  a  uniform  upper  bound  for  and  ITj^  and  therefore  the  unifomi 
exponential  stability  of  =  T^-B^F^and  S^  =  T^-F^C^.  We  conclude  that  nf.jPj.^  and 

^  y\ 

converge  strongly  in  V  to  ITj^  and  ITi^  ,  respectively  .and  that  FnPn  and  Fj^  converge  to  F 
and  F  in  norm  as  N  -^  <».  The  approximating  optimal  functional  feedback  control  and  observ'er 
gains,  fj^  and  fj.^,  converge  respectively  to  f  and  f  in  the  norm  as  N 

In  implementing  the  approximation  scheme  just  outlined  above,  eigenvector  decomposition  of 
the  associated  Hamiltonian  matrix  was  used  to  solve  the  matrix  algebraic  Riccati  equations  (3.10)  and 
(3. 1 1)  (see  Pappas,  et.  al.,  1980).  The  required  matrix  exponentials  also  were  computed  using 
eigenvalue/eigenvector  decomposition.  All  calculations  were  carried  out  via  Fortran  codes  on  an 
IBM  PC  AT.  Weseta  =  V.l,q  =  q=  r  =  r=1.0,  ^  =  V2/2  andx  =  .01  and  obtained  the 
functional  gains  plotted  in  Figs.  4.1  and  4.2 .  We  plot  fj.^  and  fj.^  as  well  as  Dfj,^  and  Dfj^  to  exhibit 
the  convergence.  We  note  that  Df  (or  Dfj.j)  appears  as  the  feedback  kernel  in  the  optimal  control 
law  (4.7). 

We  also  simulated  the  operation  of  the  closed-loop  system  with  an  approximating  compensator. 
Using  a  20  mode  model  for  the  infinite  dimensional  system  and  N  =  12,  we  computed  the  closed- 
loop  spectrum  of  the  approximating  compensator  (i.e.  the  eigenvalues  of  the  operator  given  by 

(3.9)  with  N  =  12).  These  eigenvalues  along  with  the  first  20  open-loop  eigenvalues  (i.e.  the  first 
20  eigenvalues  of  the  operator  T  =  tr(x))  and  the  approximating  closed-loop  control  and  observer 
eigenvalues  are  tabulated  in  Table  4.1  below.  Table  4.1  reveals  that  the  last  seven  open-loop 
eigenvalues  remain  essentially  unchanged  in  the  closed-loop  system-i.e.  these  modes  are  neither 
controlled  nor  obsen'ed  by  the  finite  dimensional  compensator.  Also,  as  one  would  expect,  o{/6y^) 
consists  essentially  of  the  union  of  o(Sj^) ,  0(5^^)  and  the  eigenvalues  corresponding  to  the 


& 


uncontrolled/unobsen'ed  modes  of  the  open-loop  system. 


It  is  worth  noting  thn*.  the  scheme  we  have  outlined  above  for  the  Neumann  boundary  control 
problem  is  the  same  scheme  that  one  would  ordinarily  use  if  the  problem  were  formulated  in  the  space 
H  -  i.e.  if  the  output  operator  C  was  bounded  on  [>2(0.1)  (see  Gibson  and  Rosen,  1986).  This  is 
possible  primarily  because  the  space  V  =  H^CO,!)  is  the  natural  energy  space  for  the  underlying 
homogeneous  or  open-loop  system.  Consequently,  the  inherent  self-adjointness  and  coercivity  in  the 
problem  is  preserved  when  it  is  formulated  in  the  stronger  space.  In  the  case  of  Dirichlet  boundary 
control,  the  situation  is  quite  different. 

For  the  Dirichlet  boundary  control  system  (4.1),  (4.2),  (4.4)  -  (4.6),  we  choose  the  spaces  H,  V 

and  W  and  their  corresponding  inner  products  to  be  the  same  as  they  were  in  the  Neumann  case. 

The  operators  A  6  J5(W,H)  and  C  e  33(V,R*)  also  remain  unchanged,  however  now  we  have 

r  e  S(W,R^)  given  by  Fq)  =  <p(l).  It  then  follows  that  the  operator  Ci  ;  Dom(Cl)  c  H  — >  H  is  given 

by  Ci9  =  aD\  for  cp  e  H^(0,1)  fl  (0,1).  It  is  well  known  that  Ci  is  densely  defined,  negative 

definite  and  self-adjoint  and  that  it  is  the  infinitesimal  generator  of  the  uniformly  exponentially  stable 

analytic  semigroup  { tr(t) :  t  ^  0}  of  bounded,  self-adjoint  linear  operators  on  H.  However  this 

time  the  operators  O' (t)  for  t  >  0  are  neither  self-adjoint  nor  a  semigroup  on  V.  Indeed,  since 

R(0' (t))  c  (0,1)  for  all  t  >  0  and  (0,1)  is  a  closed  proper  subspace  of  (0,1),  O' (t)  is  not 

strongly  continuous  in  the  V-norm  at  t  =  0.  (The  fact  that  our  general  framework  requires 

IT^  =  1  and  lfi.(r^)  c  V  precludes  our  choosing  V  to  be  (0,1).)  On  the  other  hand, 

0 

{ O'(t) :  t  >  0}  an  analytic  semigroup  implies  (see  Pazy,  1983)  that  there  exists  a  constant  p  >  0  foi 

t 

which  ICiO'  (t)lj^  ^  pt'^  for  t  >  0  .  Consequently,  if  we  define  T  =  O'  (x),  then  it  follows  that 
T  6  2(V)  and  moreover,  that 


IT^<^(pl^ 

V 


=  -  a-Uci0'(kx)9,0'(kx)9>  ^  a'^  I  CiO'(kx)9lH  IO'(kx)9lH 


akx 


akx 
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for  k  =  1,2,...  and  cp  e  V.  We  have  therefore 


(4.9)  IT’^Iv  =  l(T*)'^lv  Mr^  k  =  0. 1 ,2.... 


where  M  >  0  and  r  <  1 . 

We  again  choose  e  2j(R^W)  as  (  P''u)(x)  =  xu  for  x  e  [0,1].  Then  c;  TICA)  and  we 

have  reformulated  the  boundary  control  system  (4.1),  (4.2),  (4.4)  -  (4.6)  in  the  general  form  of  (2.1) 
-  (2.4)  and  conditions  1)  -  5)  are  satisfied. 

We  formulate  the  optimal  control  problem  with  the  perfomiance  index 


J(u)=  £  q<w^,w^>„  +  ru; 
where  q  >  0  and  r  >  0.  That  is,  we  take  Q  to  be  the  bounded,  self-adjoint  nonnegative  operator  on 

1  fX  ft 

Hj^(0,l)  given  by  (Q<p)(x)  =  q  J  g  J  y  <P(z)tizdy  and  R  to  be  r.  For  the  estimator  problem  we  set 

^  Aik  Aik  As 

Q  =  q  I  and  R  =  r  with  q  >  0  and  f  >  0. 

The  uniform  exponential  bound  (4.9 )  implies  the  existence  of  unique,  nonnegative,  self-adjoint 

As 

solutions  n  and  FI  to  the  algebraic  Riccati  equations  (3.1)  and  (3.5).  The  optimal  control  is  again  of 

^  1 

the  form  (4.7)  with  the  optimal  functional  gains  f  and  f  in  . 

The  fact  that  { «r(t) :  t  >  0}  is  not  a  semigroup  on  V  precludes  the  use  of  a  semigroup  -  theoretic 
approach  to  approximation.  We  therefore  employ  modal  subspaces  and  approximate  the  open-loop 
state  transition  operator  T  directly  as  a  bounded  linear  operator  on  V. 

For  each  N  =  1,2,...  let  Vj^  =  span  where  for  x  e  [0,1],  cpoW  “  ^  9iW  =  sinjnx, 

j  =  1 ,2,...,N.  Let  Pf.j  denote  the  orthogonal  projection  of  H  =  L2(0, 1 )  onto  span  { cp; } 

J  j  *■  1 

and  let  Pj^  denote  the  orthogonal  projection  of  V  onto  Vn-  Using  the  fact  that  V  =  H^(0,1)  ©  (pQ,  it 
is  not  difficult  to  see  that  Pj^cp  =  cp(l)9o  +  Pn^9  "  9(1)9o)  for  cp  e  V  and  hence,  via  elementary 
properties  of  Fourier  series  (see  Tolstov,  1962),  that  i(Pf^  -  Reply  =  l(pj^  -  I)((p  -  (p(l)(po)|y  ->  0 
as  N  oo  for  each  cp  e  V. 

We  define  T^  e  S(Vj^)  by  T^  =  P^T.  Then,  since  ?l(T)  =  R(0r(x))  c  (0,1) , 


for  Vn  =  S  ‘Pj  ^ 
j=o 

Vn  =  Pn'^n  =  =  Pn^ (^)Vn  =  ^ ('^)PnVn  =  S  {  +  Vn)  e' ''  >j . 

j=i 

It  follows  that  T^  =  Pn^*’  *V  ~  k  =  0,1,2,...  with  M  >  0  and  r  <  1  independent 

of  N,  and  that 


KTnPn  '  T)(plY  ^  l(Pj,jTPj^  -  Pj^T)(ply  +  l(Pj^  -  I)T(pIy 

<  Mrl(Pj^  -  I)(pIy  +  I(Pn  -  I)T(plv  ^  0 


as  N  — >  oo  for  tp  e  V.  Similarly,  T'  strongly  on  V  as  N  ->  <». 

The  approximating  input,  output,  and  state  penalization  operators  Bj.^,  ^N*  Qn  and  Qn  take  the 


B^u  =  (I  -  Tj^)r^u  =  V  +  X  —  "Pj'^  ‘ 

j=l 

=<1PnQ  Qn~  Reasoning  as  we  did  in  the  Neumann  case,  the  approximating 
algebraic  Riccati  equations  (3.7)  and  (3.8)  admit  unique,  nonnegative,  self-adjoint  solutions  and 

^  /N  ^  ^ 

ITjy  respectively,  ITf^Pj*^  —>  IT  and  Hf^PN  —>  IT  strongly  on  V  and  F^Pn  — >  F  and  Fj^  — >  F  in  norm 
as  N  ^  oo .  The  approximating  functional  feedback  control  and  observer  gains  fj^  and  f^^  converge 

A 

to  f  and  f  respectively,  strongly  in  H*  as  N  . 

With  a  =  1.0,  q  =  q  =  r  =  1.0,  f  =  5.0,  ^  =  "il/l  and  x  =  .01  and  the  scheme  outlined  above  we 
obtained  the  approximating  optimal  functional  feedback  control  and  observer  gains  plotted  in  Figs. 
4.3  and  4.4  below.  The  first  12  open-loop  and  the  approximating  closed-loop  control  and  observer 
eigenvalues  for  N  =  12  are  tabulated  in  Table  4.2. 

Table  4.2  reveals  an  interlacing  of  the  closed-loop  control  and  open-loop  eigenvalues.  That  is, 
the  closed-loop  control  eigenvalues  (i.e.  the  elements  in  the  spectrum  of  S)  are  alternately  more  and 
less  stable  than  the  corresponding  open-loop  eigenvalues.  We  also  have  observ'ed  this  phenomenon 
in  other  numerical  studies  we  are  carrying  out  involving  LQG  boundary  control  for  flexible 
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structures.  In  additon,  in  the  Dirichlet  boundary  control  system  discussed  above,  if  Q  is  cliosen  as 
the  identity  operator  on  V  =  virtually  all  of  the  closed-loop  control  eigenvalues  are  less 

stable  than  the  corresponding  open-loop  eigenvalues.  It  is  clear  that  this  non-standard  behavior 
results  from  the  presence  of  the  one  dimensional  subspace  represented  by  irKP*").  Indeed,  the 
behavior  of  the  closed-loop  spectrum  in  the  case  of  Neumann  boundary  control  is  as  would  be 
expected.  We  feel  that  what  we  are  seeing  can  most  likely  be  explained  via  infinite  dimensional 
analogs  of  existing  results  relating  the  asymptotic  properties  of  the  closed-loop  spectrum  of  a  linear 
regulator  and  the  zeros  of  the  corresponding  open-loop  transfer  function  (see  Kwakemaak  and  Sivan, 
1972  and  Harvey  and  Stein,  1978).  However,  as  of  yet,  we  have  been  unable  to  establish  this 
conjecture  satisfactorily  and  we  consider  it  to  be  beyond  the  scope  of  this  paper,  which  is  primarily 
concerned  with  approximation.  We  leave  it  as  an  interesting  open  question. 


5.  Concluding  Remr 


We  have  developed  a  framework  for  the  finite  dimensional  approximation  of  optimal  discrete-time 
LQG  compensators  for  distributed  parameter  systems  with  boundary  input  and  unbounded 
measurement.  Our  theory  applies  to  the  class  of  boundary  control  problems  which  can  be  formulated 
in  a  state  space  in  wliich  both  the  discrete-time  input  and  output  operators  are  continuous.  We  have 
used  a  functional  analytic  treatment  to  develop  a  convergence  theory  and  have  demonstrated  the 
feasibility  of  our  approach  via  examples  involving  either  the  Neumann  or  Dirichlet  boundary  control  of 
a  one  dimensional  heat  equation  with  point  measurement  of  temperature.  We  have  shown  that  while 
both  problems  outwardly  appear  to  be  quite  similar,  they  in  fact  require  very  different  approaches  to 
approximation.  Also  in  the  Dirichlet  case  the  observed  behavior  of  the  resulting  closed-loop  spectrum 
is,  in  some  ways  unexpected  and  its  explanation  remains  open. 

Finally,  we  have  been  looking  at  the  application  of  our  schemes  to  LQG  problems  for  flexible 
structures  with  boundary  inputs  and  unbounded  measurement  and  systems  with  control  and/or 
observations  delays.  We  have  been  considering  vibration  suppression  for  cantilevered  beams  via 
shear  or  moment  inputs  at  the  free  end  and  pointwise  observation  of  strain  or  acceleration.  These 
studies  are  currently  underway  with  the  results  to  be  ref)orted  elsewhere. 
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ABSTRACT 


A  numerical  approximation  scheme  for  the  esumation  of  functional  parameters  in 
EuIer-BemouUi  models  for  the  ta.nsverse  vibration  of  flexible  beams  with  tip  bodies  is  developed. 
The  method  permits  the  identification  of  spatially  varying  flexural  stiffness  and  Voigt- Kelvin 
viscoelastic  damping  coefficients  which  appear  in  the  hybrid  system  of  ordinary  and  partial 
differential  equations  and  boundary  conditions  describing  the  dynamics  of  such  stmctures.  An 
inverse  problem  is  formulated  as  a  least  squares  fit  to  data  subject  to  constraints  in  the  form  of  a 
vector  system  of  abstract  first  order  evolution  equations.  Spline-based  finite  element 
approximations  are  used  to  finite  dimensionalize  the  problem.  Theoretical  convergence  results  are 
given  and  numerical  studies  carried  out  on  both  conventional  (serial)  and  vector  computers  are 
discussed. 
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We  develop  here  numerical  approximation  methods  for  the  estimation  of  functional  or  more 
precisely,  spatially  varying  parameters  that  describe  material  propenies  in  continuum  models  for 
elastic  structures.  In  particular,  we  consider  tlie  identification  of  the  flexural  stiffness  and 
Voigt-Kelvin  viscoelastic  damping  coefficients  in  Euler-Bernoulli  models  for  the  transverse 
vibration  of  long,  slender,  flexible  beams  with  tip  appendages.  The  primary  motivation  for  the 
study  we  report  on  here  is  the  modeling  and  ultimately  the  control  of  the  dynamics  of  large  flexible 
spacecraft.  The  type  of  structures  to  which  we  are  referring  includes  satellites  with  flexible 
appendages  (solar  panels  and  the  like)  antennas  (reflectors  as  well  as  supporting  structures)  and 
trussed  masts  and  platforms,  both  shuttle  attached  and  free  flying. 

The  difficulties  involved  in  the  design  of  efficient  and  practical  control  laws  and  in  particular  the 
need  for  extremely  high  fidelity  models  for  stmetures  of  these  types  are  well  documented  (see,  for 
example,  [1],  [8],  [21],  [22]).  Their  high  flexibility,  light  damping,  construction  wtith  new  and 
relatively  untested  composite  materials  (usually  graphite-epoxy)  and  overall  complexity  together 
with  their  use  in  a  fuel  limited  and  highly  variable  environment  all  contribute  to  making  space 
structure  stabilization  and  control  a  formidable  task.  It  is  becoming  increasingly  clear  that  the  use 
of  continuum  or  distributed  models  with  spatially  and  /  or  temporally  varying  functional  parameters 
has  the  potential  to  offer  several  distinct  and  significant  advantages.  Included  among  them  is  the 
ability  to,  in  some  sense,  capture  the  physics  and  inherent  infinite  dimensionality  of  the  dynamics 
while  at  the  same  time  greatly  reducing  the  number  of  unknown  or  experimentally  indeterminable 
material  parameters  which  have  to  be  identified  (see  [15],  [18],  [23],  [28],  [35]). 

In  our  study  we  have  considered  exclusively  Voigt-Kelvin  viscoelastic  damping  which  is  based 
on  the  hypothesis  that  the  damping  moment  is  proportional  to  strain  rate.  There  exists  considerable 
evidence  to  suggest  that  damping  mechanisms  in  composite  materials  are  significantly  more 
complex  than  the  one  described  by  the  Voigt-Kelvin  model.  For  example,  it  has  been  conjectured 
by  some  investigators  that  an  appropriate  model  might  involve  hysteretic  or  hereditary  effects. 
However,  since  there  are  a  number  of  materials  for  which  the  Voigt-Kelvin  assumption  is 
appropriate  and  moreover,  since  at  present  many  questions  regarding  the  modeling  of  structural 


damping  mechanisms  remain  open,  we  feel  that  the  Voigt-Kclvin  model  leads  to  a  reasonable  class 
of  examples  and  problems  on  which  we  can  begin  to  develop,  test,  and  evaluate  identification 
schemes. 

Our  treatment  here  is  similar  in  spirit  to  some  of  our  earlier  efforts  and  the  work  of  others  on 
inverse  problems  for  elastic  structures  (see  [2],  [3],  [4], [5],  [6],  [14],  [17],  [26],  [31]). 

Formulating  the  identification  problem  as  a  least  squares  fit  to  data,  the  scheme  we  develop 

I 

involves  a  spline  based  finite  element  approximation  to  die  hybrid  system  of  coupled  ordinary  ar.d  ' 

partial  differendal  equations  describing  the  dynamics  of  the  structure  together  with  a  spline  based 
discretization  of  the  admissible  parameter  set. 

I 

Our  approach  here  specifically  differs  from  the  one  taken  in  [5],  [6]  in  that  the  present  schem.e  is 
derived  from  an  alternative  state  space  formulation  for  the  underlying  dynamical  equations.  We  ] 

consider  the  higher  order  analog  of  the  classical  conservative  formulation  for  a  second  order 

I 

hyperbolic  equation  as  a  first  order  vector  system  in  the  natural  states  of  strain  u^  and  velocity  Uj^.  ! 

We  have  considered  identification  schemes  based  upon  this  formulation  previously  in  [31].  i 

However  by  replacing  the  semigroup  theoretic  convergence  arguments  used  there  with  weak  or 

I 

variational  arguments  (in  the  spirit  of  those  commonly  found  in  the  finite  element  literature)  as  used 
in  [5],  we  are  able  to  significantly  weaken  the  hypotheses  necessary  to  ensure  convergence.  W'e 

point  out  below  that  the  weakening  of  these  hypotheses  has  both  theoretical  and  computational  ; 

I 

significance. 

Along  with  reporting  theoretical  convergence  results,  we  discuss  numerical  findings.  Our 
computational  results  are  based  upon  extensive  numerical  studies  which  involved  a  variety  of 
examples  and  two  machines.  In  addition  to  testing  our  scheme  on  a  conventional  serial  computer 
(an  IBM  3081)  we  vectorized  our  codes  for  the  Cray  1-S  and  then  benchmarked  some  of  our  runs 
in  order  to  explore  the  potential  of  vector  architectures  in  the  context  of  inverse  problems  for 
systems  described  by  distributed  parameter  models. 

We  provide  a  brief  outline  and  summary  of  the  remainder  of  the  paper.  In  Section  2  we  specify 
the  ordinary  and  partial  differential  equations  which  govern  the  underlying  dynamics  of  the 
structure  and  precisely  formulate  the  identification  problem.  We  reformulate  the  initial-boundary 
value  problem  as  an  abstract  second  order  evolution  equation  and  then  as  a  first  order  vector 
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system.  Hxistence,  viniqucncss  and  regularity  results  for  solutions  arc  summarized.  Section  3 
contains  the  abstract  approximation  theory  and  convergence  results.  A  spline-based  scheme  is 
discussed  in  detail  in  Section  4  and  our  numerical  findings  are  reported  and  summarized  in  Section 
5. 

We  use  standard  notation  throughout.  For  X  and  Y  Banach  spaces,  tlie  Banach  space  of 
continuous  linear  transfonnations  from  X  into  Y  is  denoted  by  jS  (X,Y).  When  X  =  Y  we  use  the 
shortha.  .lOtation  S  (X).  The  spaces  of  (equivalence  classes  oO  functions  f  from  an  interval  d 
into  X  which  satisfy 

J  I  f(6)  I  ^  d9  <  OO  or  ess  sup  !f(,0)|^<°o 
d  ^ 

are  denoted  respectively  by  L2(d  ;  X)  and  L,^(d  ;  X).  For  k  ==  0,1,2...  the  space  of  X-valued 
functions  with  k  continuous  derivatives  on  d  are  denoted  by  C’^(d  ;  X).  When  k  =  0  we  use 
C(d  ;  X).  The  completion  of  the  space  C’'(d  ;  X)  with  respect  to  the  norm 

ifb  =  (^zp’(9)q'd0^ 

is  denoted  by  H^(d;  X),  When  X  =  R  we  use  simply  L^Cd),  L^(d),  C’^Cd)  and  H’'(d). 

2.  The  rdentification  Problem  ^ 

We  consider  the  identification,  or  estimation,  of  the  mass  and/or  material  properties  of  a  long, 
slender,  flexible,  viscoelastic  beam  of  length  Z  and  spatially  varying  mass  density  p  which  is 
clamped  at  one  end  and  free  at  the  other  with  a  body  rigidly  attached  at  the  free  end  (see  Figure  2. 1 


Figure  2. 1 


We  assume  chat  the  material  behavior  of  the  beam  is  that  of  an  idealized  Voigt-Kelvin  solid  wiuh 
modulus  of  elasticity  E  and  coefficient  of  viscosity  Cq  (see  [30  ]).  We  assume  further  that  E, 
and  the  cross  sectional  moment  of  inertia  I  of  the  beam  are  in  general  spatially  varying.  We  take  the 
mass  properties  of  the  tip  body  to  be  mass  m  and  moment  of  inertia  J  about  the  center  of  mass  O 
which  is  assumed  to  be  located  at  a  distance  c  from  the  tip  of  the  beam  directed  along  the  beam's  tip 
tangent  (see  Figure  2.2  below). 


Figure  2.2 
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We  note  that  there  is  no  essential  loss  of  generality  in  assuming  that  the  mass  center  of  tlie  tip  body 
is  not  offset  from  the  tip  tangent  of  the  beam.  We  refer  the  interested  reader  to  [31]  where  ilie 
more  general  situation  is  treated.  Also,  the  problem  with  non-zero  mass  center  offset  can  be 
transformed  into  a  problem  of  the  general  fonn  of  the  one  which  will  be  considered  here.  Sec  [32] 
for  details. 

Letting  u  =  u(t,x)  denote  the  transverse  displacement  of  the  beam  at  position  x  at  time  t  and 

assuming  only  small  deformations  ( |  u(t,x)  |  «  -6,  |  —  (t,x)  |  «  1),  the  Euler-BemouIIi  theory 

dx 

a:.  !  elementary  Newtonian  mechanics  yield  the  hybrid  system  of  ordinary  and  partial  differential 
equations  (see  [19],  [34]) 
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3  u  a  3  u  3  u 

(2. 1)  p^t.x)  (El— (t,x)  -f-  CoI-^(t,x) }  = 

3x  3x  3t 


3 

~  o—^bx)  +  f(t,x), 
3x  3x 


xe(0,-6),  t  >  0 


(2.2)  iii^(t,i  )  me  (tj)  -  (EI^  -I- 

9t  3t  3x  3x  3x 

3^u 

Cq! — r-)(b-^)  =  o-^t,-i)  +  g(t),  t  >  0 

3x^3t 

2 

(2.3)  rnc^-^tj)  -h  (J  -h  mc^)  (tJ)  -h  EI-^  (tJ)  -h 

dt  3t  3x  3x 


^  3^u  du 

CqI — (t,i)  -  -ca^tj)  +  h(t), 

3x  3t 


(2.4)  u(t,0)  =  0,  — (t,0)  =  0.  t>0 

3x 


(2.5)  u(O.x)  -  <;)(x). 


(O.x)  =  v(x). 


X£[0.^]. 


Equation  (2.2)  and  (2.3)  arc  derived  from  the  usual  transverse  and  rotational  equilibrium 
considerations  at  the  free  end.  The  geometric  boundary  conditions  (2.4)  are  the  zero  displacement 
and  zero  slope  constraints  at  the  clamped  end.  The  functions  f  =  f(t,x),  g  =  g(t),  h  h(t)  and 
a  =  a(t,x)  denote  externally  applied  loads  in  the  fonn  of  moments  (h)  and  transversally  (f  and  g)  or 
a.\iaUy  (a)  directed  forces  exerted  on  the  beam  or  tip  body.  (In  fact,  h(t)  =  h(t)  +  cg(t)  where  h  is 
an  externally  applied  torque  on  the  tip  body).  The  temporal  boundary  conditions  (2.5)  reflect  the 
initial  displacement  and  veloaty  distributions  which  are  assumed  to  be  given  by  the  functions  (J)  and 


V/  respectively. 


We  treat  the  initial-boundary  value  problem  (2.1)  -  (2.5)  in  the  form  of  an  abstract  second  order 


evolution  equation  which  we  then  rewrite  as  an  equivalent  tirst  order  vector  system.  The  particular 
state  space  formulation  we  choose  forms  the  basis  for  the  finite  dimensional  approximation 
schemes  we  develop  in  the  next  section.  It  also  allows  us  to  easily  establish  existence,  uniqueness 
and  regularity  results  for  solutions  to  (2.1)  -  (2.5)  using  the  theory  of  abstract  parabolic  systems. 

Let  H  denote  the  Hilbert  space  x  1^(0, 1)  with  inner  product 
(Tl2>^2-92>H  =^lTl2  +  ^1^2  +  <9l.02>O 


and  let  V  denote  the  Hilbert  space 

V  =  { (ri.^,0)  e  H :  0  e  h2(0,  ^),  0(0)  =  D0(O)  =  0,  n  =  0(  e ),  ^  =  D0(  e ) } 
with  inner  product 


<§l,  §2>v  =  <EI(D20i),  d202>o 


for  -  (0j(-6 ),  D0j  (■£),  0j)  £  V,  i  =  1,2.  In  the  above  definitions  the  inner  product  <-,  >o  is 


the  standard  one  on  LriO,  I )  and  D  denotes  the  spatial  differentiation  operators  —  or  With  H  as 


the  pivot  space,  we  obtain  the  usual  dense  embeddings  V  c  H  =  H'  c:  V  . 


We  consider  the  system  (2.1)  -  (2.5)  in  the  form  of  the  abstract  second  order  initial  value 


IIKWM 


(2.6)  mo  Uu(t)  +  GoU^(t)  +  OCouCO  =  Bo(t)u(t)  +  ^oCO.  t  >  0 

(2.7)  G(0)  =  iT^(O)  =  V 


in  the  state  u(t)  =  (u(t,-&),Du(t,-8.),u(t,  •))  e  H.  The  abstract  mass,  damping  and  stiffness 

operators  mQ.  ^  given  formally  by 

mo(Tl,^.9)  =  (mri  +  mc^,  mcrj  +  (J  +  mc^)q,  p0) 

=  (-D(CDr(D20))(^).  CDI(D20)(^),D2(CDI(D20))) 

and 

0Co(Tl,^.e)  =  (-D(EI(D20))(^).  EI(D20)(i),  d2(EI(D“0))) 
respectively.  For  each  c>  0,  the  operator  valued  function  ?3q  and  input  or  forcing  function  D^q 
take  on  the  values 

Bo(t)(Ti,q.0)  =  (-a(t,^)(D0(i)),  -  co(t,^)(D0(e)),  D(a(t,-)(D0))) 

and 

yoW  =  (g(t).h(t).f(v)). 

The  initial  conditions  $  and  ^  are  given  by 

^=(<pa).  D<p(Z),  (}>) 

and 

Y=(Y(^),  D\|r(^),\j/). 

The  formal  definitions  given  above  can  be  made  precise  and  the  existence  and  uniqueness  of 
solutions  to  the  initial  value  problem  (2.6),  (2.7)  can  be  established  if  we  make  the  following 
assumptions. 

Aj  The  functions  p,  El  and  C^I  are  elements  in  C[0,  •£]  and  there  e.xists  a  positive  constant  a  for 
which  p(x)  >  a,  EI(x)  >  a,  C^Kx)  >  a,  x  e  (0,  -t]. 


A2  The  mapping  t  o(t,-)  is  an  element  in  L^((0,T);  H’(0,-£))  for  some  T  >  0. 

A3  The  mapping  t  is  an  element  in  ^((O.T):  ^(O.-t))  and  g,h  £  L2(0,T). 

A4  The  function  (J)  is  an  element  in  H^(0,-&)  with  (J){0)  =  D(t)(0)  =  0  and  \/  e  L2(0,-6)  with  v(^) 
and  Dy(-&)  defined. 

Under  the  hypotheses  Aj  -  A^  above,  the  operator  TIq  is  a  bounded  linear  operator  from  H  onto  H 
andtjQ:  Dom  (Uq)  c  H  — »  H  and  OCq:  Dom  (OCq)  c:  H  — >  H  are  densely  defined,  nonnegative, 
self-adjoint  operators  defined  on  Dom(GQ)  =  {9  e  V  :  CjjKD^G)  e  H^(0,-&)}  and  Dom  (0<Iq)  =  {9 
e  V:  EI(D^9)  e  H“(0,-t)}  respecdvely  (see  [32]).  For  each  t  £  (0,T),  ^  2^(V,H)  and  £ 

H  while  $  E  V  and  y  £  H.  It  also  follows  that  ?3q  £  Lj,^((0,T);  S(V,H))  and  D^q  £  L2((0,T);  H). 
We  shall  call  a  mapping  t  — >  u(t)  from  [0,T]  into  H  a  strong  solution  to  (2.6),  (2.7)  if 

u  £  C  ([0,T];  V)  n  Cl((0,T];  V)  n  Ck[0,Tl;  H)  n  C2((0,T];  H), 

u(t)  £  Dom  (OCq),  U[(t)  £  Dom  (^jq),  t  e  (0,T1,  and  u  satisfies  (2.6)  and  (2.7)  where  the  time 
derivatives  are  interpreted  in  a  strong  (norm)  sense  in  H.  We  shall  call  a  mapping  t  u(t)  from 
[0,T1  into  H  a  weaJc  solution  to  (2.6),  (2.7)  if 

Q  £  C([0,T];  V)  n  Hl((0,T);  V)  n  Ck[0,T];  H)  n  h2((0,T);  V) 

and  it  satisfies  the  initial  value  problem  (2.6),  (2.7)  with  the  operators  Uq  and  Xq  replaced  by  their 

natural  extensions  to  operators  in  3j(V,V')  and  the  time  derivatives  are  interpreted  in  a  weak  or 
distributional  sense  (see  [20  ],  [27]).  A  function  u  =  u(t,x)  will  be  called  a  strong  (weak)  solution 
to  the  initial-boundary  value  problem  (2.1)  -  (2.5)  if  the  mapping  t  — >  u(t)  given  by 
u(t)  =  (u(t,-i),Du(t,i),  u(t,-))  is  a  strong  (weak)  solution  to  (2.6),  (2.7). 

Our  approximation  theory  for  the  estimation  problem  to  be  developed  below  is  based  upon  the 

reformulation  of  the  initial  value  problem  (2.6),  (2.7)  as  a  first  order  vector  system.  This 
reformulation  is  formally  equivalent  to  rewriting  the  initial-boundary  value  problem  (2.1)  -  (2.5) 


u» j||«  iL  k 9^  * 


as  a  first  order  system  in  the  states  D^ii  (strain)  and  u,  (velocity)  (see  [3],  [31]).  We  note  that 
since  the  stiffness  operator  Xq  is  nonnegative  and  self-adjoint  it  has  a  unique  nonnegative, 

m 

selfadjoint  square  root  %q  :  V  c  H  ->  H.  It  can  be  written  in  factored  form  as 

OCo  =  LeiL 

where  L  :  V  c  H  is  given  by 

L9  =  D^e, 

for  §  =  (e(^).D0(i),  0)  e  V.  and  L|i :  Dorn  (L^j)  c  ^(0,^)  ^  H  by 


Dorn (l4)=  (9 e  1^2(0. :  EI0eH^(O,e)} 

(2.8)  •  L|i0  =  (-D(EI0)(^).  EI0(e),  d2(EI0)). 

If,  for  t  e  C[0,-ft]  with  t(x)  >  a  >  0,  x  e  [0,-E],  we  let  L2  denote  the  Hilbert  space  L2(0,'£) 
endowed  with  the  inner  product 

<01,02>o^^  =  <X0j,02>o 

then  given  by  (2.8)  with  El  replaced  by  x  is  the  Hilbert  space  adjoint  of  L  as  a  mapping  from 
V  c  H  into  L2  . 

We  note  that  L  e  S(V,L2  £j)  is  a  Hilbert  space  isomorphism  with 


and  L’^  :  1^(0.  •^)  V  given  by 


We  also  have 


•*  *  i  .X 

’^9  =  (JJ  0(y)dydx,  J0(x)dx,  J  J  0(y)dydx). 


^0  “  • 

»  L> 


Letting  “Mi  =  1^(0, -t)  x  H  with  inner  product 

(2.9)  <(8i.(i)i  ,^i.Xi)).(e2.(ri2.^2-X2))>W  =  +  <TIIo(t1i.^i.Xi).  (Tn2>^2-X2)>n 

and  V  =  1-2(0, ■&)  x  V  with  inner  product 

<(9i.Xi).(02'X2)^‘U‘  ^01-02^0, El  ^Xl.X2^V 
we  have  the  dense  imbeddings  V  c  %  c  V’.  We  consider  the  initial  value  problem  for 

z(t)  =  (w(t),  v(t))  eW  given  by 

(2.10)  w^(t)  =  Lv(t) 

(2.11)  Trig  C  (t)  =  -  w(t)  -  L^^jL(;(t)  +  (t)L'V'(t)  +  ^^(t)  0  <  t  <  T 

(2.12)  w(0)  =  L$.  C(0)  =  v 

which  we  rewrite  as 

(2.13)  z^(t)  =  G(t)z(t)  +  m  0<t^T. 

(2.14)  z(0)  =  Zo 
where 

(2.15)  a(t)  =  Cl  4- 73  (t) 

with  Cl  :  i9  c  W  — >%,  Tie  L,„((0,T);  £(H)),  IX  e  L2((0,T);  %)  and  zq  6  IHi  given  by 

^(0,X)  -  (Lx.  -  TT^o^L^e  -  TUq  L*^jLx) 
for  (9,x)  e  J9  =  Dom(L£j)x  Domf^Q), 

^f5(t)(0,(Tl,q.x))  =  (O.-Rq  73Q(t)L'’0), 

for  (9,(Tl,q,x))e  W, 

IX (t)  =  (O.Ro  IX^d)) 
and 

Zq  =  (L^.v) . 


In  formulating  the  inverse  problem,  we  keep  technical  details  to  a  minimum  by  considering  on 


the  estimation  of  the  beam's  spatially  varying  flexural  stiffness  El  and  viscous  damping  coefficient 
CqI.  Extending  the  finite  dimensional  approximation  methods  and  corresponding  convergence 
tlieory  which  arc  developed  below  so  as  to  be  applicable  to  the  identification  of  other  structural  or 
input  parameters,  for  example  mass  properties  (of  the  beam  and/or  tip  body),  initial  conditions  or 
loading,  is,  at  least  in  principle,  routine  (see  (7]  [  12]  [14]  [16  ]  [31]). 

Let  9)  =  C[0,'6]  X  C[0,(1]  with  norm 

(2.16)  Iqlq,  =  I(qi.q2)^  =  ki  loo+  IqiIoo 

=  sup  lq^(.x)l  +  sup  lq2(x)l. 

n£[e.t]  x£[0,f] 

We  taJee  the  admissible  parameter  space  Q  to  be  a  compact  subset  of  (compact  with  respect  to  the 
metric  topology  induced  by  the  norai  (2.16)).  Recalling  assumption  (i)  we  assume  further  that  the 
set  Q  has  the  property  that  all  q  =  (qj,q2)  e  Q  sadsfy  qi(x)  >  a  and  q2(x)  >  a  ,  x  e  [O,-?.]. 

We  formulate  the  identification  problem  as  a  least-squares  fit-to-data  over  the  admissible  parameter 
space  Q.  We  assume  that  the  structure  has  undergone  a  time  varying  elastic  deformation  in 
response  to  the  initial  conditions  described  by  ^  and  v/  and  the  input  loads  represented  by  f,g,h  and 
a.  Denoting  the  observation  space  by  Z,  we  assume  that  at  times  t^,  i  =  l,2,...,v  measurements 
(t^)  e  Z  (e.g.  displacement,  velocity,  slope,  strain,  etc.)  were  taken  from  the  structure. 

We  require  that  Z  be  a  linear  space  endowed  with  a  norm  I  •  I  2[  ^  denote  an 

appropriately  defined  continuous  mapping  from  W  into  Z  .  For  example,  suppose  that 
displacement  measurements  have  been  taken  at  the  points  xj,  j  =  l,2,...,|i  along  the  span  of  the 
beam.  We  choose  Z  as  Euclidean  }j.-space,  Rd »  and  take  F  to  be 

r(z)  =  (0(xjX  e(.X2),...,9(x^))’^ 

where  z  =  (w,v)  e  %  and 

0  =  (e(i’).  De(e).e)  =  L-^wev. 

With  distributed  strain  or  velocity  observations,  we  would  take  r(z)  =  w  or  r(z)  =  v  respectively. 
We  formulate  the  identification  problem  as  follows 


1 J 


(ID)  Given  zX,\  =  Find  q*  £  Q  which  minimizes 

m  =  ^  Ir(z(t.;q))-C(t)i2 

i=l 

where  for  each  q  =  (qi.q2)  £  Qi  z  (•  :  q)  =  (w(-  ;  q),  v(-  ;  q))  is  the  solution  to  the  initial  value 
problem  (2.13),  (2.14)  or  (2.10)  -  (2.12)  with  El  set  equal  to  q^,  and  C^I  set  equal  to  (\2- 

It  is  immediately  clear  that  the  optimization  problem  given  above  is  inherently  infinite 
dimensional.  The  admissible  parameter  set  Q  is  a  subset  of  a  function  space  and  the  evaluation  (and 
therefore  minimization)  of  the  least-squares  performance  index  ;]  requires  the  solution  of  an  infinite 
dimensional  evolution  equation.  The  introduction  of  finite  dimensional  approximations  is  essential 
to  the  development  of  practical  computational  methods.  Fundamental  to  the  approach  we  take  here 
is  a  weak,  distributional,  or  variational  formulation  of  the  initial  value  problem  (2. 13),  (2. 14).  We 
derive  the  weak  form  and  briefly  outline  existence,  uniqueness  and  regularity  results  for  solutions. 

In  the  usual  manner,  we  extend  the  operator  Cl(t)  given  by  (2. 15)  to  an  operator  in  Z  (Tf,'!/') 
via 

(Cl(t)(v))(v)  =  a  (t)(v,v),  v,v  e  IT 
where  the  bilinear  form  a  (t)(-,-) :  V  x  V  ^  R  is  given  by 

(2.17) 

X 

ca(t,i)  J  0^(x)dx(Dx2('^))  “  J  ^^0 ' 

0  0 

Standard  estimates  can  be  used  to  demonstrate  the  existence  of  positive  constants  k,  X  and  P  for 
which 

la(0(vi,V2)l  ^  k  Iv^lf^;  lv2lv  ,  Vj  e  TT,  i  =  1,2, 
and 

a(t)(v,v)-hX|v|5^>p|v|^,  veV,  tF[0,T]. 

Consequently  (see  [27  ])  the  system  (2.13),  (2.14)  interpreted  as  an  initial  value  problem  in  1/'  or 
equivalently,  written  in  weak  form  as 


(2.18)  <Zi(t),v>^  =  a(t)(z(t).v)  +  <  y(0.v>5<; 


V  e  V,  t  e  (0,i] 


(2.19)  z(0)  =  zo 

admiis  a  unique  solution  z  witli  z(t)  £  V,  t  £  (0,T]  and 

z  £  L2((0.T);  V)  n  C([0.T];  %)  n  H’((0.T):  V'). 

If  z  =  (w,v)  is  the  unique  solution  to  (2.18),  (2.19)  then 

ir(t)  =  L-'w(t).  t£[0.T] 
is  a  weak  solution  to  (2.6),  (2.7)  and  it  is  unique. 

Under  somewhat  stronger  hypotheses  than  those  given  in  and  Aj  above,  the  e.xistence  of 
strong  solutions  can  be  established.  Indeed,  if  in  addition  to  A^  and  A^,  we  assume 

Aj  The  mapping  t  o(t,-)  is  an  element  in  C*([0,T];  H*(0,-6))  for  some  T  >  0 

A3  The  mapping  t  —>  f(t,-)  is  an  element  in  C*([0,T];  1^(0, -fc))  and  g,h  E  C'(0,T]  (in  fact. 
Holder  continuity  will  suffice,  see  [29],  [37]) 

then  the  family  of  operators  (Cl(t)}te[0,T]  generates  a  unique  evolution  system 

(U(t,s)  ;  0  ^  s  ^  t  T)  on  IHl  and  z  given  by 

t 

(2.20)  z(t)  =  U(t,0)Zg  +  J  U(t,s)Cr(s)ds,  0  <  t  ^  T, 

0 

is  the  unique  solution  to  the  inital  value  problem  (2. 13),  (2. 14)  and  satisfies  z(t)  £  t  £  (0,T] 
with  z  £  C([0,T];  %)  n  Ck(0,T];  W) .  Once  again,  with  z  =  (w,v)  now  given  by  (2.20), 
u(t)  =  L‘'w(t),  t  £  [0,T],  is  a  strong  solution  to  (2.6),  (2.7)  and  it  is  unique. 


•A  •««  yt 


3.  An  Abstract  Approximation  Framework 

We  turn  next  to  a  discussion  of  a  general  approximation  framework  and  convergence  tlieory  for 
the  identification  problem  (ID)  fomiulated  above.  In  the  following  section  we  formulate  a  specific 
spline-based  scheme  to  which  the  general  theory  developed  here  applies. 

Fundamental  to  our  approach  is  the  construction  of  a  sequence  of  finite  dimensional  (with  regard 
to  both  the  state  dynamics  and  the  admissible  parameter  set)  approximating  identification  problems 
each  of  which,  under  appropriate  hypotheses,  can  be  shown  to  have  a  solution  that  in  some  sense 
(specifically,  subsequential  convergence)  approximates  a  solution  to  the  original  infinite 
dimensional  estimation  problem  (ID). 

In  the  discussion  to  follow,  we  exhibit  the  explicit  dependence  on  q  =  (qi,q2)  e  of  the  % 
inner  product  <  and  the  bilinear  form  a(t)(-,-)  given  in  (2.9)  and  (2.17)  respectively  by 
using  the  notation  <-,->q  and  a(t;q)(-,-)-  For  each  =  1,2,...  and  each  N2  =  1,2,...  let 

W  and  V  be  finite  dimensional  subspaces  of  L^(0,.^)  and  V  respectively.  If,  for 

N  Ni  >^1  N 

N  =  (N^,Np  we  define  V  =  W  x  V  ,  then  *1/  is  a  finite  dimensional  subspace  of  both 
W  and  V  .  Let  P  ^ V  ^  denote  the  projection  map  of  %  onto  V  ^  given  by 

(3.1)  F  (w,v)  =  (p;'w,  PjV) 

N  ^1  V 

where  P^  is  the  orthogonal  projection  of  L2(0,i)  onto  W  and  is  the  orthogonal 

projection  of  H  onto  V  ,  both  computed  with  respect  to  the  standard  (unweighted)  inner 
products  on  the  respective  spaces  1^(0, i)  and  H. 

The  Galerkin  equations  in  TT^  corresponding  to  the  system  (2.18),  (2.19)  and  q  e<^  are 

(3.2)  <z[^(t),  v^>q  =  <3(t ;  q)(z^(t),v'^)  -1-  <!J'(t),  .  v^  e  V  ,  0  <  t  <  T 

(3.3)  z‘’^(0)  =  (P%  . 


For  each  Mj  e  =  {l,2,...},i=  1,2,  let  Sj^^^be  finite  dimensional  subspaces  of 

C[0,-&]  and  for  M  =  (Mj,  M2)  define  c  by  x  Let  and 

1  2 

denote  mappings  from  C[0,/?]  onto  and  respectively  and  define  dj^,  a  mapping  from 
<1  onto  by 

I  2 

^M(q)  =  q  =  q2)  £  ^  • 

We  define  a  sequence  of  approximating  admissible  parameter  spaces  {Qj„i}.  Me  x  Z.^.  by 

(3.4)  Qm  =  ^  m(Q) 

and  formulate  the  sequence  of  approximating  identification  problems  as  follows: 

(ID,^j)  Given  ^(t.)  e  Z,i  =  l,2,...,v,  find  (q^^)  e  which  minimizes 


/(<0  =  2)  lr(z''(t|;q))-«t)|2 


i=l 


over  Qj^,  where  z^(-;  q)  is  the  solution  to  the  initial  value  problem  (3.2),  (3.3)  in 


N  -N  ^  i  *-M  •  ^ 

We  choose  bases  {9.  }•_,  ,  (x-  }-_t  f  and  {\/  }.J,  for  the  finite  dimensional 

1  1“  i  I  I”  I  IVX  1-^  I  i 


N  N2  ^ 

*  f  r  ^  o  * 


1  _ 


spaces  W  ,  V  ,  and  Sj^  respectively.  Then  e  ^  solution 


z^(  • ;  q^^)  to  the  initial  value  problem  (3.2),  (3.3)  with  q  =  qj^  =  (qj,^,j ,  qj^ )  can  be  written  as 


••vv>  < 


t  e  [O.T], 


N 

z 


(t;qM)  = 


i=l  i=l 


respectively.  Moreover,  Z^(- ;  is  the  solution  to  the  initial  value  problem  in 
R  given  by 

(3.5)  m^'^(aj^)Z^(t)  =  A'^Ct ;  Z^(t)  +  P\t).  t  e  (O.T] 

(3.6)  Z^(0)  =  z5^. 

Here  the  positive  definite  matrix  'ITl^(0Cj,^)  is  of  the  form 


m^(V 


Tn.^(aM)  = 


o 


where  7n*^(aj^)  is  a  K^-square  matrix  with  components 

N  V  t  M  N 

[T^lCVlij= 

k=l 

and  Tn.2  is  K2 -square  matrix  with  entries 

[m^iij  =  <moxr.£y>H. 

For  each  I  >  0  Ihe  matrix  A^(.t ;  a^)  is  given  by  A% ;  a^)  =  A^(CIm)  with 

i 
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E''(aM) 


^^«m)  = 


M  T  NT 

-eV,)‘  -cV,) 


B^t)  = 


D^(t) 


where  E^(aj^)  is  a  Kj^  x  matrix  with  components 


tE''(V)ii=E«i,<4>M9r.D"xr>„. 


C^Cccm)  is  a  -square  matrix  with  components 


[C  =  2L  Oj^  <Vj^D  Xi  .  D  Xj 

k=l 


and  D^(t)  is  a  x  matrix  with  components 


[D^(t)]jj  =  ca(t.i)  J  0‘'^(x)dxDxf(i)  -  <a(t ,  •)  J  0j^(x)dx.  Dxf>Q  . 

0  0 

The  nonhomogeneous  term  F^(t)  is  given  by  F^(t)  =  (0,  F^(t))  where  F^(t)  is  a  vector  with 
entries  * 


-N  , 


[F^(t)]i  =  <0^0(0, 

The  initial  data  is  of  the  form 


z^  =  (gV(4.4)' 


MM  MM 

where  the  IC^  vector  and  the  vector  are  given  component-wise  by 


[z-oiii = er>„ 


,  -  -N 


respectively,  and 


with  G[^  a  K|^- square  matrix  defined  by 


[Gflij  =  <eN  ef>o . 


and  G^  a  K^-square  matrix  with  components 


=  <x\^ ,  xf  >H  • 


ir 

It  is  now  easily  seen  that  the  finite  dimensional  identification  problem  (fDj^)  in  fact 

involves  simply  the  minimization  of  a  least-squares  performance  index  over  a  subset  of 

N 

R  .  Furthermore,  the  evaluation  of  the  functional  ;)  requires  only  the  solution 

N  N 

of  the  Kj  -I-  dimensional,  linear,  noa-autonomous  ordinary  differential  equation  (3.5) 


wilh  initial  conditions  specified  in  (3.6).  If  the  existence  of  solutions  to  the  finite  dimensional 

optimization  problems  can  be  established,  it  is  immediately  clear  that  they  can,  in  principle,  be 

computed  using  standard  techniques.  Conditions  which  guarantee  the  existence  of  solutions  to 
N 

problem  and  the  fact  that  they  in  some  sense  approximate  solutions  to  the  original 
infinite  dimensional  estimation  problem  (ID)  are  given  in  the  followng  theorem. 

Theorem  3.1  Suppose 

Hj  the  mappings  are  continuous  from  Q  into 

H2  for  each  q  e  Q,  — >  q  as  I  M  I  — >  with  the  convergence  being  uniform  in  q  for  q  e  Q, 

H3  the  spaces  and  projections  '3P^  are  such  that  if  (q"’^}  is  a  sequence  in  Q  with 

qM  ->  q  =  (q^,  q2)  e  Q  as  I  N  I  «  then  z^(t ;  q^^)  ->  z(t ;  ^  in  L2(0,-6)  x  H  for  each  t  e 
[0,T]  as  I N  1  «  where  z^(-  ;  q^)  is  the  solution  to  the  initial  value  problem  (3.2),  (3.3) 

with  q  =  q^  and  z(- ;  q)  is  the  solution  to  the  initial  value  problem  (2.18),  (2.19) 
corresponding  to  El  =  qj^  and  C^I  =  q2. 

Then,  each  of  the  problems  (IDj_^)  has  a  solution  (q,^)  .  Furthermore,  the  sequence  { (qj^p  } 
admits  a  <^-convergent  subsequence  whose  limit  q  is  a  point  in  Q  and  is  a  solution  to  problem 

(ID). 

In  the  statement  of  the  theorem,  for  an  element  K  =  (K^,K2)  e  x  we  have  adopted  the 
notation  I K  I  ->  00  to  denote  K^.K^  <»  .  We  remark  that  it  is  also  true  that  the  limit  point  of 

any  ^(-convergent  subsequence  {(q  ^  1  of  1  with  I M  1 ,  |  N  |  -><»  as  j3c  ^  is  a 

solution  to  problem  (ID)  as  well.  Moreover,  if  problem  (ID)  has  a  unique  solution,  q*,  then  the 

N  •  • 

sequence  {(qj^)  }  itself  converges  to  q  .  It  is  also  important  to  note  that  the  hypotheses  of  the 
theorem  do  not  require  that  c  Q. 

We  have  established  results  analogous  to  those  given  in  Theorem  3.1  for  inverse  problems 
involving  parabolic  and  hyperbolic  systems  (see,  for  example  [12],  [13],  [16])  as  well  as  for 
related  methods  for  higher  order  equations  for  elastic  structures  (see  [4],  [5],  [6]).  For  the  flexible 
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structure  problems  treated  here,  the  essential  features  of  the  argument  remain,  for  the  most  part, 
unchanged.  We  therefore  only  briefly  sketch  them  below. 

Standard  continuous  dependence  results  for  linear  ordinary  differential  systems,  the  continuity 
assumptions  on  and  T  (and  therefore  on  as  well)  and  the  fact  that  Q  is  a  compact  subset  of 
H  are  sufficent  to  conclude  that  there  exists  a  solution  e  to  problem  (ID^^  . 

-M 

The  definition  of  the  space  (see  (3.4))  implies  the  existence  of  a  e  Q  for  which 

N  •  — N  — M 

(Qm)  “  Q  compact,  there  exists  a  subsequence  { q  .}  of  { q^^j }  with 

-N^  •  .  -N^ 

q  •  — >  q  eQasjJc^oo.  The  subsequence  {q  .}  can  always  be  chosen  with 
1  mJ  1 ,  1 1  — >  oo  as  j,k  — >  oo  .  It  follows  that 


d  ((q  j))^d  (q). 

M 


and  consequently  that 


n’' 

d  ((q  i))^d  (d  :(q)). 


qeQ. 


Assumption  H2  above  and 


(q  ;)  — >  q  as  j Jc  00  .  Taking  the  limit  as  jjc  — > in  (3.7)  above  with  an 

application  of  assumpdori  H3,  we  find  !l(q  )  ^  <I(q),  q  £  Q  ,  and  hence  that  q  is  a  solution  to 
problem  (ID). 


In  this  section  we  outline  a  scheme  which  uses  piecewise  polynomial  spline  functions  and  show 
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f  ^  j  j  f 


that  it  satisfies  tlie  conditions  and  hypotheses  of  Theorem  3.1.  We  first  treat  the  discretization  of 
the  admissible  parameter  set  Q. 

I  2 

For  each  M  =  (Mp  .Ny  eZ^xZ^  let  and  Aj^.j  denote  the  uniform  partitions  of  the  interval 
[0,-&]  determined  by  the  meshes  {0,-6/Mp  2-E/Mj, . . .  ,1]  and  2-6/M2,. 

respectively.  For  m  =  1,2,.  . .  and  A  a  partition  of  [0,  i]  let  Sp(m,A)  denote  the  usual  spline 

space  of  functions  in  which  are  polynomials  of  degree  2m- 1  on  each  subinterval  of 

i  i  i 

A  (see  [36]).  We  then  define  =  Sp(l,  A),  i  =  1,2.  In  this  case  we  have  dim  =  Mj  + 

i 

1,  i  =  1,2,  with  the  usual  "hat"  functions  fonning  a  cardinal  basis  for  each  of  the  spaces  S  ,  i  = 

i  i  ’ 

1,2.  For  i  =  1,2,  let  C[0,i]  S  ,^,  be  the  interpolation  operator  defined  by 

(v)(^)=v(4)’  j  =  oa,2 . M, 

M.  i 

for  Y  e  C[0,-6].  The  theory  of  interpolatory  splines  (see  [33])  yields  the  continuous  dependence 
result 

1^mYi-^Ly2Ioo^  lYl-Y2loo,i=1.2 

where  Yp  Yj  £  C[0,lt]  and  consequendy  that  hypothesis  Hj  of  Theorem  3.1  is  satisfied.  Also, 
the  approximation  result  (see  [36]) 

N’mY-yL  ^  w(Y,  1/M.) 

where  a)(Y,5)  is  the  usual  modulus  of  continuity  of  ye  C[0,Z]  with  respect  to  5,  together  with  the 
assumption  that  Q  is  a  compact  subset  of  =  C[0,-6]  x  C[0,-^]  and  die  Arzela-Ascoli  theorem 
yield  that  hypothesis  is  satisfied  as  weU. 

Next  we  define  a  state  approximation  and  verify  that  hypothesis  holds.  As  above,  given 
N  =  (Np  N2)  e  X  ,  we  define  the  uniform  partitions  Af^of  the  interval  [0,i]  determined  by 
the  meshes  [0,-i/N-,  2i/N-,.  i  =  1,2.  We  may  then  choose  either 


In  the  first  case,  once  again  the  "hat "  functions  may  be  chosen  as  a  basis  with 
^1  N 

dim  W  =  Kj  =  Nj  +  1.  In  the  second,  the  standard  cubic  B-splines  (see  [33]), 

Ni  Nj^-l  ^ 

{Bj  }j  _  1  ,  corresponding  to  the  partition  fonn  an  appropriate  finite  element  basis  with 

dim  W  =  +  3.  In  either  case,  approximation  results  for  interpolatory  splines  can  be  used 

to  obtain 

(4.1)  IPf^e-elg^OasN^ -^oo 

for  0  6  L2(0,-t). 

We  set 


=  ((X(i),  D%(i),x)  e  H  :  X  e  Sp(2,A^^ ),  x(0)  =  Dx(0)  =  0}  . 

^2 

Then  V  c  V  and  defining 


N  N  N  N 

2  ^2  ^2  ^2 


Pi  -2B/.2B/. 

^2  ^2 

P.  =  ,  i  =  2,3, . .  .  :N2  +  1 


and 


N. 


K. 


p.  ^  =  (P.  ^(i),Dp.  ^(i).p.  \  i  =  1.2 . +  1. 


N, 


N, 


N. 


the  collection  (P. forms  a  basis  for  V  ^  with  dim  V  ^  =  K^  =  N2+  1.  With  V  ^ 
as  defined  above,  it  is  not  difficult  to  show  (using  arguments  similar  to  those  in  [31]) 


(4.2)  Ip?(t1.^,X)-(t],^,X)Ih-^0  asN2 
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for  (ri.^.X)  e  H  and 


LP^X-Lx1h->0  asN2-^oo 


for  xeV- 


So  as  to  avoid  obscuring  the  essential  features  of  our  argument  with  technical  details,  we  verify 
hypothesis  H3  for  the  spline-based  scheme  described  above  in  the  case  a  h  0.  The  term  which 
results  from  axial  loading  is  a  bounded  perturbation  and  does  not  involve  the  unknown  parameters. 
Showing  that  the  desired  convergence  continues  to  hold  in  the  presence  of  a  non-zero  axiady 
directed  acceleration  requires  only  a  routine  extension  of  the  proof  which  we  give  below  (see  [14]). 

Suppose  that  {q^}  is  a  sequence  in  Q  with  q*^— >qeQ  as  |n|  — ><».  Let  z^  =  (w^, 
denote  the  solution  to  (3.2),  (3.3)  with  q  =  q^  and  let  z  =  (w,'^  denote  the  solution  to  (2.18), 
(2.19)  corresponding  to  q.  We  shall  require  the  assumption  that  z  is  a  strong  solution. 

In  the  estimates  which  follow,  we  simplify  our  notation  by  referring  to  the  inner  products 
(norms)  <•,•>  ^(M  )  and  <•,•>- (|- 1  -)  on  W  by  <-,->^  ( |  •  I  )  and  <•,•> 

q  q  IN 

( I  •  I )  respectively.  Also  note  that  with  CT  =  0,  we  have  a  (t;  q)(-,-)  =  a  (q)(-,  ). 


Since 


llz^  -  zll  ^  llz^  -  yNzIl  -I-  ll(F^  -  I)zll 


where  n  il  denotes  the  usual  (unweighted)  product  norm  on  %=  L2(0,'E-)  x  H,  (3.1),  (4.1)  and 
(4.2)  imply  that  we  need  only  to  consider  the  term  llz^  -  '3’^zll .  Letting  y‘’^(t)  =  z^(t)  -  '3P^z(t), 
using  (2. 1 8),  (2.19),  (3.2),  (3.3)  and  the  fact  that  c  V  we  find 


N  N  /  ^  N  .  N  N 

>N  =  <(Z  -  ^  z)j ,  V  -I-  <z^ ,  V  >  -  <Zj ,  V  >J^ 


(3(q^)(y^,v^)  -  a  (q^)(z  -  ?’^z,v^)  +  a(q^)(z,v^)  -  (3(q)(z,v^) 
+  <0f',v%  -  <CF',vN>  e  VN  0  <  t  ^  T 

yN(0)  =  0. 


Choosing  e  we  obtain 


=  <(z  -  iP^z)^ ,  +  <(qj  -  q,^)Wt ,  w^  -  P,^w>o 

-  <q|^L(v  -  P2  v),  -  P,^w>o  +  <q{^(w  -  P|^w),  L(\^  -  P^ v)>g  - 

<q^L(v  -  P^v),  L('/^  -  P2v)>o  +  <(q|^  -  qi)Lv,  w^'^-  Pj^w>o  - 
<(qi^  -  qi)w,  L(v^  -  P^^>0  -  <(q^  -  q2)Lv,L('?^  -  P^^>o  - 

Recalling  that  Q  is  a  compact  subset  of  q,  and  that  for  q  =  (q^,  q2)  e  Q.  El  =  qi  and  CqI  =  q2  axe 
assumed  to  satisfy  assumption  Aj  of  Section  2,  we  find 

1  lly^ll^  +  1  L(C''^-  p^v)  i  Ko{ll(I .  ip\/ 

+  lly^ll^+  Iqj-q^lMwJj  +  (w^-P^wl2  + 

lLa-P^)vlo+ Iw^-Pfwlj  +^|(I-Pf)w|j 

+  e  iu^  -  p^;)lj  +  ^  |La-P'2)v|j  + 
elL(v^-P^v)|J+  lqJ^-qJ^lLvlJ  + 

Bl^v"-F^C)|J.^|q^qJilLv|5 

+  elL(v'^-P2;)l2} 

where  iqj  is  a  positive  constant  and  e  is  an  arbitrary  positive  constant  Gathering  up  like  terms 
and  choosing  e  <  “  we  obtain 


Ai|yV<K-^{ll(I-/)z/+  I\vj2+  lL(r-P2)vl2  + 


Ia-P^)wl5  +  lq^qlll  1Lv|J+  iq^qililwl2+  |q^  -  qj2  1  L^l^}  +  K^HyV 

where  Kj  and  ^2  are  positive  constants.  Integrating  both  sides  of  the  above  inequality  from  0  to  t 
and  recalling  (4.4)  we  obtain 


(4.5)  |ly^(t)ll^<5  +  icj||y%ll^ds 


where 

T 

5  =  <1 1 (li(I  -  )z^(s)ll^  +  I  ^  1  i  I  I J  +  1  L(I  -  P2)v(s)  I J  +  I  (I  -  P^)w(s)  I J  + 

0 

1^2  ■  ^2  111  • 


Since  ->  q  as  1 N 1  ^  and  z  =  (w,v)^  was  assumed  to  be  a  strong  solution,  (4.1),  (4.2), 
(4.3)  and  (4.5)  together  with  an  application  of  the  Gronwall  inequality  yield  the  desired  result, 
l!y‘'^(t)II->0. 

A  close  inspection  of  the  estimates  above  reveals  that  they  depend,  to  a  large  extent,  on  the 
presence  of  the  viscous  damping  term  <CqI  Lx^,  LX2>o  bilinear  form  a  (t)(-,-)  given  in 
(2. 17).  That  is,  we  require  that  q2  ^  a  >  0  for  some  a  >  0  for  all  q  =  (qp  q2)  e  Q.  In  the 
absence  of  the  Voigt-Kelvin  damping  we  can  still  argue  the  convergence  of  z^  to  z;  however,  we 
must  assume  that  Q  is  tf-compacL  If  one  is  to  enforce  the  compactness  constraint  on  Q  when 
solving  the  finite  dimensional  optimization  problems  (a  desirable  implementation  feature  in  many 
cases  -  see  [10],[1 1]),  this  stronger  assumption  becomes  especially  unappealing.  On  the  other 
hand,  by  employing  a  somewhat  different  (but  closely  related)  factorization  of  the  stiffness 
operator  Xq  than  the  one  which  was  used  here  (one  which  is  formally  equivalent  to  rewriting  the 
initial- boundary  value  problem  (2.1)  -  (2.5)  as  a  first  order  system  in  the  states  El  D^u  and  Uj  as 
opposed  to  and  Uj)  hypothesis  H3  of  Theorem  3.1  can  be  verified  for  the  resulting  spline-based 


scheme  under  the  present  assumptions  on  Q.  Unfortunately  this  scheme  is  also  difficult  to 
implement  and  from  a  numerical  standpoint,  has  not  perfonned  as  satisfactorily  as  the  one  based  on 
the  formulation  given  in  this  paper.  The  present  scheme  performed  well  whether  or  not  damping 
was  present  in  the  equation  and  hence  the  assumption  that  C^I  >  a  >  0  may  be  an  artifact  of  our 
proof  of  convergence  (see  Example  5.3  below). 

5.  Numerical  Findings 

We  present  and  discuss  some  of  the  results  which  we  obtained  from  our  numerical  studies  of  the 

scheme  that  was  described  in  Section  4.  All  codes  were  written  in  FORTR_\N,  and  tested  and  run 

on  the  EBM  3081  at  either  Brown  University  or  the  University  of  Southern  California.  The  same 

codes  were,  with  only  minor  modification,  run  on  the  Cray  1-S  at  Boeing  Computer  Services  in 

Seattle  with  support  made  available  to  us  through  the  National  Science  Foundation's  Super 

Computer  Initiative  program.  Examples  were  benchmarked  so  that  the  potential  benefits  of 

vectorization  to  our  research  program  could  be  accurately  and  effectively  assessed.  Our  findings 

are  described  below.  This  information  will  become  especially  important  to  us  when  we  begin  to 

consider  the  extension  of  our  general  approach  to  inverse  problems  involving  the  vibration  of  two 

dimensional  structures,  such  as  flexible  plates  or  platforms,  or  vibrations  of  structures  in  which 

nonlinearities  play  a  significant  role.  The  finite  dimensional  optimization 
N 

problems  (DDj^)  were  solved  using  the  IMSL  routine  ZXSSQ,  an  implementation  of  the  iterative 

Levenberg-Marquardt  quasi-Newton  algorithm.  The  finite  dimensional  initial  value  problems 
(3.5),  (3.6)  were  solved  in  each  iteration  of  the  minimization  procedure  (for  the  evaluation  of  the 
least-squares  performance  index  iJ  and  its  gradient  v-ith  respect  to  the  parameters)  using  Gear's 
method  for  stiff  systems  (IMSL  routine  DGEAR). 

Our  codes  were  written  to  take  full  advantage  of  the  banded  structure  of  the  generalized  mass, 
stiffness  and  damping  matrices  afforded  by  the  use  of  polynomial  B-spline  elements.  All  necessary 
inner  products  were  computed  using  a  two  point  composite  Gauss- Legendre  quadrature  scheme. 

All  of  the  examples  presented  here  involve  fits  based  upon  displacement  measurements  obtained 
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llirough  simulation.  "True"  values  (which,  in  the  examples  below  will  be  denoted  with  an  asterisk, 
for  example  El  ,  CqI  ,  etc.)  for  the  unknown  parameters  were  chosen.  The  resulting  initial 
boundary  value  problem  (2.1)  -  (2.5)  was  then  solved  using  an  independent  integration  scheme. 
(\Ve  used  a  seven  element,  quintic  spline  based  Galerkin  method  applied  directly  to  the  second 
order  system  (2.6),  (2.7)).  This  procedure  produced  sufficient  noise  in  the  data  so  that  the  use  of  a 
random  noise  generator  was  not  required. 

In  addition  to  the  test  example  numerical  studies  we  report  on  here  we  have  successfully  used 
methods  similar  to  those  developed  above  with  experimental  data.  These  results  are  presented  in 
detail  in  [9]. 

In  the  examples  which  follow  we  took  the  axial  loading  to  be  induced  by  an  acceleration  of  the 
base  or  root  of  the  structure  in  the  positive  .x-direction.  In  this  case  we  have  (see  [34]) 


a(t,x)  =  -3^(0  [m  +  J  p(y)dy} 


where  m  is  the  mass  of  the  tip  body,  p  is  the  linear  mass  density  of  the  beam  and  a^  is  the  time 
dependent  base  acceleration. 

In  Examples  5.1  thru  5.4  below  we  took  I  =  \,  p(x)  =  3-  xfor  O^xSl,  f(t,x)  =  e^sin  27tt, 
g(t)  =  2e*^  h(t)  =  e’2^  aQ(t)  =  1  for  0  S  t  <  1.5,  aQ(t)  =  0  for  t  >  1.5,  m  =  1.5,  c  =  .1  and  J  =  .52 
and  considered  the  estimation  of  the  flexural  stiffness  coefficient  El  and/or  the  viscoelastic  damping 
coefficient  CqI  only.  In  Example  5.1, 5.2  and  5.4,  the  fits  we  describe  are  based  upon 
observations  at  times  t^  =  .2i,  i  =  1,2, ...  ,5  at  locations  xj  =  .5,  .75  and  1.  In  Example  5.3 
observations  at  times  t^  =  .5i,  i  =  1,2,...,10  at  locations  xj  =  .75  and  1  were  used.  In  all  of  the 

N 

examples  we  discuss  here  the  space  W  was  generated  by  cubic  splines  (i.e.  as  Sp  {2,A^ )) 

with  =  N2  =  N.  This  corresponds  to  the  approximation  of  the  first  and  second  components  of 

z  with  respectively  N  +  3  and  N  +  1  piecewise  cubic  elements. 

The  compacmess  constraints  on  the  spaces  Qm  were  not  enforced  when  the  finite  dimensional 

N 

optimization  problems  (IDj^)  were  solved.  When  Mj  and  M2  became  large,  the  inherent 
ill-posedness  of  the  inverse  problem  became  apparent  as  the  performance  of  our  schemes 

i 

deteriorated.  There  is  evidence  strongly  suggesting  that  this  situation  can  be  remedied  either  by 


IMJIUMM  U  Jl  JfWiJULV  mj  J.m  AV.V.V.  J.  j.  K 


'In 


imposing  the  compactness  constraints  on  the  admissible  parameter  space  and  then  solving  the 
minimization  problem  using  a  constrained  optimization  procedure  (see  [10],  [1  Ij)  or  by 
regularizing  the  least  squares  performance  index  (see  [24],  [25]).  We  intend  to  direct  our  attention 
to  these  ideas  in  the  near  future. 

Example  5.1 

In  this  example  we  consider  the  simultaneous  estimation  of  a  constant  flexural  stiffness 
coefficient,  El*  =  .15,  and  a  damping  coefficient  given  by  CqI*(x)  =7(1.5  -  tanh  (3x  -  1.5)), 
xe  [0,1],  with  7=  .01.  With  N  =4,  Mj  =  I  and  M2  =  3  and  taking  start  up  values  (for  the  least 
squares  minimization  algorithm)  EI^  =  .1  and  CqI°(x)  =  .015,  x  e  [0,1]  we  obtained  the  results 
shown  in  Figure  5.1  below.  This  particular  run  required  approximately  30  seconds  of  CPU  time 
on  the  IBM  3081 


Figure  5. 1 

We  observed  that  how  weU  the  scheme  performed  depended  upon  the  magnitude  of  the  scaling 
factor  y.  As  y  was  decreased,  so  too  did  the  "sensitivity"  of  the  least  squres  performance  index  to 
the  damping  coefficient  Results  similar  to  those  shown  in  the  figures  above  were  obtained  with 
y  =  .005.  With  y=  .001,  on  the  other  hand,  we  were  unable  to  simultaneously  identify  both  of  the 

unknown  parameters.  However,  again  with  y=  .001,  but  this  time  fixing  El  at  the  true  value,  we 
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were  able  to  identify  CqI  alone. 

*  *  j 

When  we  replaced  the  constant  El  with  the  linear  function  El  (x)  =  1  -  —  x  and  took 
Y=  1,  the  performance  of  the  scheme,  from  a  qualitative  point  of  view,  remained  unchanged. 

Example  5.2 

We  again  consider  the  simultaneous  estimation  of  the  stiffness  and  damping  coefficients.  We 
again  set  El*  =  .15  but  this  time  choose  CqI*(x)  =  .01  (1.5  -  tanh  (20x  -  10)),  x  e  [0,1].  The 
identification  of  this  steeper  hyperbolic  tangent  function  has,  in  past  test  examples,  proven  to  be  a 
somewhat  stiffer  challange  for  our  methods  (see  [5], [6]).  With  N  =  4,  =  1,  M2  =  3,  EI*^  =  .1 

and  CDl°(x)  =  .015for0<x<  1,  we  obtained  the  estimates  which  are  plotted  along  with  the  true 
parameters  in  Figure  5.2. 


Figure  5.2 

Also,  although  the  theory  was  not  explicitly  treated  here,  we  note  that  elements  other  than  linear 
splines  can  be  used  to  discretize  the  admissible  parameter  space.  Our  investigations  have  included 
numerical  studies  with  0-order  splines  (i.e.  piecewise  constant  functions)  and  cubic  spline 
elements.  Using  two  linear  elements  to  approximate  El  (i.e.  Mj  =  1)  and  nine  cubic  elements  to 
discretize  C^I  we  obtained  the  estimates.shown  in  Figure  5.3.  We  have  obtained  an  acceptable 
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Figure  5.3 


Example  5.3 

In  this  example  we  identify  only  the  spatially  varying  flexural  stiffness  coefficient 
El  (x)  =  1.5  -  tanh  (3x  - 1.5),  x  e  [0,1],  in  a  model  with  no  viscoelastic  damping  (CqI  =  0).  In 
Section  4  we  remarked  that  our  convergence  arguments  required  either  the  presence  of  viscoelastic 
damping  in  the  model  or  that  the  admissible  parameter  set  Q  be  compact  in  the  stronger 
topology.  The  results  shown  in  the  figure  below  suggest  that  this  is  only  an  artifact  of  our  proof 
and  not  a  fundamental  requirement  for  the  convergence  of  our  approximation  (i.e.  the  absence  of 
damping  does  not  appear  to  effect  the  overall  performance  of  our  scheme). 

Taking  N  =  4  and  Mj  equal  to  1  thru  8  we  produced  the  results  shown  in  the  series  of  graphs  in 
Figure  5.4.  The  initial  estimate  or  start-up  value  for  El*  was  taken  to  be  the  constant  function 
El0(x)  =  1  forO^x^l. 

i 

Recalling  our  earlier  remarks,  the  oscillations  which  appear  in  the  graphs  corresponding  to 
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Mj  =  6,1  and  8  due  to  the  inherent  ill  posedncss  of  the  estimation  problem  are  not  unexpected. 

In  fact,  as  Mj  or  M2  — >  <»,  the  appearance  of  the  undesirable  oscillations  in  our  final  csiintatcs 
occurred  in  virtually  every  test  we  ran.  As  we  have  noted  earlier  however,  preliminary  findings  in 
related  studies  [10]  and  [1 1]  regarding  the  enforcing  of  the  compactness  constraints  and  the 
subsequent  use  of  constrained  optimization  techniques  to  solve  the  approximating  finite 
dimensional  identification  problems  suggest  that  this  difficulty  can  be  overcome.  Our 
investigations  in  these  directions  are  continuing. 

In  addition,  the  series  of  tests  corresponding  to  the  graphs  in  Figure  5.4  were  benchmarked  on 
the  IBM  3081  and  the  Cray  1-S.  The  same  estimates  were  obtained  on  both  machines.  However, 
we  were  able  to  achieve  a  speed-up  factor /,  of  7  -  10  on  the  vector  machine.  The  CPU  times  are 
reported  in  Table  5.1.  In  comparing  the  CPU  times  on  the  3081  for  this  example  with  the  times 
reported  for  the  previous  examples  it  is  important  to  note  that  the  results  here  were  based  upon 
observations  taken  over  the  longer  time  interval,  [0,5],  versus  the  interval  [0,1]  for  examples  5.1 
and  5.2. 


Example  5.4 

In  Figure  5.5  below  we  plot  the  final  estimates  obtained  when  we  attempted  to  use  our  scheme 
to  simultaneously  identify  the  spatially  varying  flexural  stiffness  coefficient 
El  (x)  =  .5  -t-  4x(l  -  x),  X  e  [0,1],  and  viscoelestic  damping  coefficient, 

Hfi 

CqI  (x)  =  .1  (1.5  -  tanh  (3x  -  1.5)),  x  e  [0,1] .  The  start-up  values  for  the  iterative  least  squares 
minimization  routine  were  taken  to  be  the  constant  functions  EI°(x)  =  1  and  CqI*^(x)  =  .15  for  0  <  x 
<  1.  The  graphs  in  the  figure  were  obtained  with  N  =  4  and  a  linear  spline  discretization  of  the 
admissible  parameter  set  Q  with  =  M2  =  3.  In  all  of  our  tests  with  this  example  the  minimum 
sum  of  the  squares  of  the  residuals  was  in  the  range  10'^  -  10'^  with  the  optimization  typically 
requiring  50  -  70  seconds  of  CPU  time. 
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Figure  5.7 

For  this  example  we  also  tried  a  cubic-spline  based  discretization  for  Q.  We  considered  all 
possible  combinations,  linear  splines  for  El  together  with  cubic  splines  for  C^I  ,  cubic  splines  for 

4t  4( 

El  together  with  linear  splines  for  CqI  ,  etc.  Although  small  values  for  the  sum  of  the  squares  of 
the  residuals  were  obtained  in  each  instance,  our  by  far  best  approximation  to  the  true  parameters  is 
the  one  shown  in  Figure  5.5  which  corresponds  to  a  linear  spline  based  discretization  for  both 
components  of  the  admissible  parameter  set. 

Holding  CqI  fixed  at  the  tme  value  and  using  cubic  splines  to  identify  El  and  then  holding  El 

fixed  at  the  tme  value  and  using  cubic  splines  to  identify  CdI*  we  were  able  to  obtain  the  estimates 
plotted  in  Figures  5.6  and  5.7  respectively.  The  estimate  for  El  graphed  in  Figure  5.6  was 
obtained  with  10  cubic  elements  while  the  estimate  for  CqI*  in  Figure  5.7  is  a  linear  combination 
of  6  cubic  elements.  An  inspection  of  these  figures  reveals  that  while  the  approximations  obtained 
are  at  least  marginally  acceptable,  it  is  also  not  surprising  that  our  scheme  had  some  difficulty  when 
we  attemped  to  identify  both  parameters  simultaneously  with  a  cubic  spline-based  discretization  for 
either  one  or  both  components  of  Q. 

For  this  example  we  also  looked  at  the  robustness  of  our  iterative  scheme  with  respect  to  the 
initial  values  chosen  (i.c.,  El®  and  Cfjl®)!  In  Figure  5.8  we  plot  those  points  in  the  Cpjl®  -  El® 


plane  which  correspond  to  the  start  up  values  we  tried.  Tlie  point  marked  with  "  *  "  corresponds  to 
the  startup  values  which  produced  the  appro.ximations  shown  in  Figure  5.5.  The  points  m.arked 


with  "  X  "  correspond  to  start-up  values  which  led  to  essentially  the  same  extimates  as  tliose  shown 

in  the  figures.  The  points  marked  with  an  "  ©  "  correspond  to  start-up  values  for  which  the 
scheme  did  not  converge.  The  region  whose  boundary  is  denoted  with  dashed  lines  corresponds  to 
a  "convergence  envelope"  for  the  vector  valued  function  (  CdI  ,  El  ).  An  analogous  study  was 
carried  out  for  Example  5.2,  for  which  similar  robustness  results  were  obtained. 


Finally  we  offer  several  summary  comments  on  some  of  our  other  numerical  findings.  In 
virtually  all  examples  we  tried,  we  found  that  the  estimates  yielded  by  the  scheme  which  we 

develop  here  based  oh  state  space  coordinates  (D^,Uj)  and  the  ones  yielded  by  the  scheme  based 

on  a  state  space  formulation  in  coordinates  (u,u^)  described  in  [5]  and  [6]  were  comparable. 
Although  in  any  given  example  one  scheme  or  the  other  may  produce  a  somewhat  better 

approximation  to  the  true  parameters,  we  found  it  impossible  to  designate  or  identify  a  clear  favorite 
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among  the  two  methods. 
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We  also  rcui  a  series  of  tests  in  which  we  varied  the  boundtury  conditions  at  the  free  end  of  the 
beam.  That  is,  in  addition  to  the  tip  body  end  condition  we  considered  a  beam  which  is  clamped  at 
one  end  and  free  at  the  other  with  either  a  point  mass  (c  =  J  =  0)  or  no  mass  (m  =  c  =  J  =  0)  rigidly 
attached  at  the  tip.  We  also  studied  the  effect  that  the  presence  or  absence  of  external  forces  and/or 
moments  at  the  tip  of  the  beam  (i.e.  g  and  h)  has  on  the  performance  of  our  scheme.  Based  upon 
these  tests,  we  found  it  difficult  to  make  definitive  statements  regarding  "best"  experimental 
procedures  for  identification  of  structural  parameters  with  our  schemes.  However,  we  are  able  to 
offer  several  observations.  For  example,  with  a  point  mass  at  the  tip,  the  schemes  performance 

was  enhanced  when  an  external  moment  was  applied  at  the  tip  (i.e.  h  ^  0).  On  the  other  hand,  the 

presence  of  an  externally  applied  force  in  the  transverse  direction  (i.e.  g  0)  did  not  appear  to  have 
any  effect  at  aU.  Also,  with  no  mass  at  the  tip,  the  scheme  was  most  effective  when 
g  =  h  =  0.  In  general  we  found  the  scheme  to  be  most  dependable  with  tip  body  end  conditions. 
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We  develop  a  computational  method  for  the  estimation  of  parame¬ 
ters  in  a  distributed  model  for  a  flexible  structure.  The  structure 
we  consider  (part  of  the  "RPL  experiment")  consists  of  a  cantilevered 


beam  with  a  thruster  and  linear  accelerometer  at  the  free  end.  The 
thruster  is  fed  by  a  pressurized  hose  whose  horizontal  motion  effects 
the  transverse  vibration  of  the  beam.  We  use  the  Euler-Bernoulli 
theory  to  model  the  vibration  of  the  beam  and  treat  the  hose-thruster 
assembly  as  a  lumped  or  point  mass-dashpot-spring  system  at  the  tip. 
Using  measurements  of  linear  acceleration  at  the  tip,  we  estimate  the 
hose  parameters  (mass,  stiffness,  damping)  and  a  Voigt-Kelvin 
viscoelastic  structural  damping  parameter  for  the  beam  using  a  least 
squares  fit  to  the  data. 

We  consider  spline  based  approximations  to  the  hybrid  (coupled 
ordinary  and  partial  differential  ecjuations)  system;  theoretical 
convergence  results  and  numerical  studies  with  both  simulation  and 
actual  experimental  data  obtained  from  the  structure  are  presented  and 
discussed. 


The  difficulties  involved  in  the  design  of  practical  and  effi¬ 
cient  control  laws  for  large  flexible  spacecraft  (e.g.  the  inherent 
infinite  dimensionalitY  of  the  SYstem,  a  large  number  of  closelY 
spaced  modal  frequencies,  high  flexibilitY*  light  damping,  a  fuel- 
limited,  hostile,  highlY  variable  environment,  etc.)  have  stimulated 
research  into  the  development  of  SYStem  identification  and  parameter 
estimation  procedures  which  will  Y^sld.  high  fidelitY  models.  A  partic¬ 
ular  area  of  interest  involves  schemes  for  the  estimation  of  material 
parameters  describing,  for  example,  mass,  inertia,  stiffness  or 
damping  properties  in  distributed  models  for  the  vibration  of 
viscoelastic  SYStems-specificallY,  mechanical  beams,  plates  and  the 
like.  In  addition,  since  the  resulting  inverse  problems  are  often 
infinite  dimensional,  substantial  attention  has  been  focused  on 
approximation;  see,  for  example,  tl],  [2],  [3],  [4],  [8]  and  [12]. 

In  these  treatments,  the  parameter  estimation  problem  is  formulated  as 
a  least  squares  fit  to  measurements  of  either  displacement  or 
velocitY-  Although  significant  gains  have  been  made  in  the  development 
of  instrumentation  to  measure  displacement  and  velocitY  (e.g.  laser 
technologY,  etc.),  one  of  the  least  expensive,  most  reliable  and  most 
commonlY  used  sensors  is  the  linear  accelerometer.  While  in  principle 
it  is  possible  to  Integrate  acceleration  measurements  once  or  twice 
to  obtain  respectivelY  velocitY  or  displacement  data,  in  practice  this 
task  can  pose  significant  challenges.  For  example,  integration  of  the 
signal  could  result  in  the  amplification  of  low  frequencY  measurement 
noise  or  dYnamic  effects  which  have  not  been  included  in  the  underlY- 
Ing  model.  In  light  of  this,  we  have  undertaken  to  show  here,  both 


theoretically  and  computationally,  that  a  scheme  in  the  spirit  of 
those  developed  in  the  previously  cited  references  can  also  be 
effectively  used  with  acceleration  measurements.  In  particular  we 
note,  this  involves  the  nontrivial  extension  of  the  familiar 
variational  arguments  which  are  used  to  demonstrate  the  convergence  of 
the  finite  element  state  approximations  upon  which  the  identification 
schemes  are  based.  Indeed,  it  must  be  shown  that  in  addition  to  the 
convergence  of  the  displacement  and  velocity,  the  convergence  of 
acceleration  can  be  obtained  as  well. 

The  other  primary  motivation  for  the  present  effort  is  that  while 
these  methods  have  been  extensively  tested  and  evaluated  with  simula¬ 
tion  data,,  they  have  never  been  tried  with  actual  experimental  data. 

We  have  tested  our  scheme  with  data  obtained  from  an  experimental 
structure  which  was  designed  and  constructed  at  the  Charles  Stark 
Draper  Laboratory  in  Cambridge,  Massachusetts  with  funding  provided  by 
the  United  States  Air  Force  Rocket  Propulsion  Laboratory  (RPL) .  The 
RPL  structure  (as  it  will  henceforth  be  referred  to  as)  was  designed 
to  serve  as  a  test  bed  for  the  implementation  and  evaluation  of 
control  algorithms  for  large  angle  slewing  of  spacecraft  with  flexible 
appendages.  The  structure  was  specifically  designed  to  exhibit 
structural  modes  and  damping  characteristics  representative  of 
realistic  large  flexible  space  structures . 

In  Section  2  we  describe  the  RPL  structure  (its  geometry, 
instrumentation,  etc.)  and  formulate  an  inverse  problem  involving  a 
distributed  system.  In  Section  3,  we  use  the  resulting  infinite 
dimensional  estimation  problem  to  motivate  the  development  of  a  finite 
dimensional,  finite  element  based  approximation  scheme.  We  also 
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discuss  our  theoretical  convergence  results.  In  Section  4  we  present 
numerical  findings. 

We  use  standard  notation  throughout.  For  X  a  normed  linear 
space.  L(X)  denotes  the  space  of  bounded  linear  operators  from  X  into 
X.  For  n  an  interval  and  k  =  0,1.2,  •••.  C  (n-.X)  denotes  the  space  of 
functions  from  n  into  X  which  are  k  times  continuously  strongly 
differentiable  on  n.  When  k  =  0  we  shall  simply  write  C(n:X).  A 
function  f  from  fi  into  X  will  be  said  to  belong  to  Lj(n;X)  if 


lnl«« 


dt  <  00.  For  k  =  0,1,2. 


,  H^(0;X)  denotes  the  completion 


of  C^(n;X)  with  respect  to  the  norm 


'j5o 


If,  in  addition,  X  is  a  Hilbert  space  with  inner  product  <  •  .  * 

H  Cn:X)is  a  Hilbert  space  with  inner  product 

=■  g^J^(t)>jj  dt. 

When  X  =•  R,  we  use  the  abbreviated  notations  C^(n),  Li(n)  and  H^(n)  . 
Note  that  H*(n)  -  Li(n)  and  is  the  standard  inner  product  on 

L,(n). 


The  RPL  structure  (see  Figure  2.1  below)  consists  of  four 
flexible  appendages  which  are  cantilevered  at  right  angles  to  one 
another  from  a  rigid  central  hub.  The  hub  is  mounted  on  an  air 
bearing  table  thus  permitting  the  near  frictionless  rotation  of  the 
structure  about  the  vertical  axis. 


FLEXI2IZ  HCSSS 


Figure  2 . 1 


Two  of  the  appendages  (which  are  mounted  to  the  hub  IBO"  apart)  are 
"active":  each  has  two  nitrogen  cold  gas  thrusters  mounted  in  opposing 
directions  at  its  tip.  The  remaining  two  appendages  are  "passive" 
with  only  counter-balancing  masses  affixed  to  their  free  ends.  The 
presence  of  the  tip  masses  on  the  passive  arms  serves  to  preserve  the 
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overall  symmetrY  of  the  structure.  Nitrogen  gas  from  tanks  mounted  on 
the  central  hub  is  supplied  to  the  thrusters  via  two  stainless  steel 
mesh-wrapped  high  pressure  hoses.  The  expulsion  of  propellant  from 
the  thruster  nozzles  is  controlled  by  electro-mechanical  or  solenoidal 
valves .  Each  of  the  four  appendages  is  equipped  with  a  sensor  in  the 
form  of  a  linear  accelerometer  attached  at  its  tip.  Data  from  the 
accelerometers  is  processed  and  recorded  and  control  input  signals  to 
the  thrusters  are  generated  by  a  MING  11/23  microcomputer.  A  detailed 
description  of  the  structure's  design  specifications  can  be  found  in 
[6]  and  [15] . 

The  problem  which  is  of  primary  concern  to  us  here  involves  the 
modeling  of  the  effects  of  the  nitrogen  supply  hoses  on  the  transverse 
vibration  of  the  active  members.  We  ccnsider  therefore,  the  structure 
with  the  central  hub  immobilized  and  look  only  at  the  vibration  of  one 
of  the  active  appendages  and  view  it  as  a  simple  cantilevered  beam 
(see  Figure  2.2). 


Figure  2.2 
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We  treat  the  thruster  assembly  as  a  point  mass  that  is  rigidly  attach¬ 
ed  to  the  beam  at  the  tip  and  propose  a  model  for  the  hose  effects  in 
the  form  of  a  proof  mass  which  reacts  against  the  tip  mass.  In 
effect,  we  consider  the  idealized,  simplified  structure  depicted  in 
Figure  2.3  below  involving  a  single,  cantilevered,  flexible,  uniform 
beam  with  a  two-mass-dashpot-spring  system  affixed  to  its  free  end. 
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Figure  2.3 


In  formulating  a  mathematical  model  for  the  structure  shown  in 
Figure  2.3  above,  we  assume  that  the  beam  is  of  length  JL  with  uniform 
rectangular  cross  section  of  height  h  and  width  b.  We  let  u(t,x)  and 
y(t)  denote  respectively  the  transverse  displacement  of  the  beam  at 
position  X  along  its  span  and  the  displacement  of  the  proof  or  hose 
mass,  each  at  time  t.  Both  are  measured  relative  to  the  x-axis  in  the 
coordinate  frame  determined  by  the  longitudinal  axis  of  the  beam  in 


its  undeformed  state  with  origin  located  at  the  beam's  root  or  fixed 
end.  Assuming  the  beam  undergoes  only  small  deformations  (i.e. 

3  VL 

|u(t,x)|  <<  t  and  |— (t,x)|  <<  1)  and  has  a  small  height  to  span 

length  ratio,  the  Euler-Bernoulli  theory  (see  [5])  including  Voigt- 
Kelvin  viscoelastic  structural  damping  (see  [10])  yields  the  partial 
differential  equation 


(2.1) 


3*u  3*  3u,  d*VL 

p  - (t.x)  +  c_I - (t.x)  +  El  - (t.x) 

3t*  ^  3x*  3t  3x* 


0. 


0  <  X  <  i , 


t 


0 


where  p  is  the  linear  mass  density  of  the  beam,  E  is  the  modulus  of 

elasticity,  c^  is  the  coefficient  of  viscosity  and  I  is  the  second 

moment  or  moment  of  inertia  of  the  cross  sectional  area  A  about  the 

neutral  axis.  For  the  beam  we  consider  here  with  constant  rectangular 

3 

cross  section,  I  =  bh  /12.  Since  the  beam  is  assumed  to  be  uniform, 
the  parameters  p,  E  and  Cj^  are  taken  to  be  constant  in  time  and  space. 

Balancing  forces  at  the  free  end,  elementary  Newtonian  mechanics 
yields  the  equations  of  motion 


(2.2) 


8‘u  3*  3u  3*u^ 

- (t,Jl)  -  c^I - (t,)l)  -  El  - (t,!L) 

^  3t*  3x*  3t  3x* 

dv  3u 

-  Cg  (— (t) - (t,i))  +  k^  (y(t)  -  u(t,i))  +  f(t), 

dt  3t 


t  >  0 


and 

(2.3) 


m 


d*y 


3u 


H 


dt  = 


-(t)  +  c„  ( — (t) - (t,!))  +  k^  (y(t)  -  u(t,i))  -  0, 

*  “  dt  ^ 


3t 


t  >  0 


for  the  tip  and  hose  masses  m^  and  mjj  respectively-  Here  is  the 
hose  stiffness,  Cjj  is  the  hose  damping  coefficient  and  f(t)  is  the 
externally  applied  force  at  time  t  due  to  the  firing  of  the  thrusters 
mounted  at  the  tip. 

Making  the  assumption  that  the  rotatory  inertia  of  the  proof  mass 
system  is  negligible,  rotational  equilibrium  at  the  tip  can  be 
expressed  as 

3*  3u  ,  a*u, 

C2.4)  C^I - (t,*.)  +  El  - (t.i)  -  0,  t  >  0. 

^  8x*  3t  ax* 

The  zero  displacement  and  zero  slope  constraints  at  the  fixed  end  are 
given  by 

3u 

(2.5)  u(t.O)  -  0  and  — (t.O)  -  0,  t  >  0 

ax 


respectively.  Taking  the  structure  to  be  initially  at  rest  we  have 
the  initial  conditions 


(2.6)  u(0,x)  -  0 
and 

(2.7)  y(0)  -  0 


— (0,x)  -  0, 

at 


— (0)  -  0. 
dt 


0  s  X  5  i 


In  the  mathematical  model  given  by  (2.1)  -  (2.7)  above  the  parameters 
p,  m,j  and  I  can  be  measured  or  computed  directly.  The  modulus  of 
elasticity  E  is  typically  determined  in  the  laboratory.  For  the  most 
commonly  used  materials  (including  aluminum  which  is  the  material  from 
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which  the  structure  of  interest  to  us  here  is  made)  its  value  can  be 
readily  looked  up  in  standard  engineering  tables.  The  parameters  Cj^, 
mjj,  Cjj  and  kjj  on  the  other  hand,  must  be  determined  experimentally; 
that  is ,  they  will  have  to  be  identified  based  upon  the  observed 
response  of  the  structure  to  a  given  input  disturbance.  This  is  one 
class  of  inverse  problems  which  we  formulate  and  consider  below.  In 
the  system  of  equations  (2.1)  -  (2.7)  we  explicitly  modeled  (albeit, 
in  a  rather  simple  fashion)  the  dynamical  effects  of  the  hose.  The 
unknown  hose  parameters  are  then  determined  as  the  solution  to  an 
inverse  problem. 

An  alternative  approach  to  obtaining  a  model  which  exhibits  a 
reasonable  degree  of  fidelity  involves  a  technique  which  is  sometimes 
referred  to  as  model  adjustment.  Starting  with  a  simple  model,  the 
parameters  are  then  "adjusted"  so  as  to  compensate  for  unmodeled 
dynamics.  The  choice  of  parameters  to  be  adjusted  and  the  resulting 
variations  may  or  may  not  be  motivated  by  physical  considerations . 

In  our  problem  for  example,  we  might  consider  a  simple  cantile¬ 
vered  beam  with  tip  mass  (i.e.  m^^  =  Cg  -  kjj  -  0)  and  then  adjust  the 
theoretical  or  measured  values  of  E  and  m^  to  compensate  for  the 
dynamical  effects  which  result  from  the  hose  mass  and  motion.  A  value 
for  the  parameter  c^  could  also  be  identified  if  damping  effects  are 
considered  significant.  Model  adjustment  was  used  in  [6]  to  obtain  a 
model  for  the  RPL  structure  upon  which  control  design  could  be  based. 

We  define  an  Inverse  problem  which  encompasses  both  of  the 
general  approaches  which  have  been  outlined  above .  We  assume  that  an 
input  disturbance  described  by  the  function  f(t),  t  €  [0,T]  is  applied 
to  the  structure  via  the  tip  thrusters  and  that  the  linear  accelera- 


tlon  at  the  free  end  of  the  beam,  z(t),  is  measured  and  recorded  for 
each  t  G  [t,,ti]  where  0  i  t,  ^  tj  i  T.  (Of  course,  in  actual 
practice,  z  could  in  fact  only  be  sampled  discretely).  Let  denote 
the  positive  real  numbers  and  let  Q  be  a  closed  and  bounded  subset  of 
R°.  We  seek  a  q  G  Q  which  minimizes 

ft  1  I  3  *u  I  2 

J(q)  -  I  - (t.iLiq)  -  z(t)  dt 

where  u(-,-;q)  denotes  the  solution  to  the  initial-boundary  value 
problem  (2.1)  -  (2.7)  corresponding  to  q  -  (m^ . E , c^ , , Cg , kg)  G  Q. 

Our  primary  concerns  in  the  next  section  will  include  well- 
posedness  of  the  system  (2.1)  -  (2.7),  existence  of  a  minimizer  for  J, 
and  development  of  approximation  techniques  to  find  this  minimizer. 


3. 
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A  computational  method  for  the  solution  of  the  estimation  problem 
posed  above  will  invariably  involve  finite  dimensional  approximation 
of  the  initial-boundary  value  problem  (2.1)  -  (2.7).  We  have  been 
successful  in  solving  inverse  problems  for  distributed  parameter 
models  for  flexible  structures  (see.  for  example,  [1],  [2],  [3],  [4], 
[12])  using  spline-based  Ritz-Galerkin  techniques.  We  apply  those 
ideas  here  and  derive  finite  element  approximations  based  upon  an 
abstract  Hilbert  space  formulation  of  the  hybrid  system  of  ordinary 
and  partial  differential  equations  and  boundary  conditions  given  in 
(2.1)  -  (2.7).  This  abstract  formulation  is  also  useful  in 
establishing  existence,  uniqueness  and  necessary  regularity  results 
for  solutions.  We  briefly  outline  the  essential  features  of  our 
general  approach  (including  theoretical  convergence  results)  in  the 
context  of  the  particular  problem  of  interest  to  us  here. 

Let  H  -  R*  X  L,(0,i)  be  endowed  with  the  usual  product  space 
inner  product 


<(^, h. 4)), (X.U, ’«')>„  -  U  +  Tiu  + 


and  let 


V  -  {(?;,n.0)  €  H:  4»  €  H*(0,JL),  (t)(0)  -  D4>(0)  -  0.  n  -  4>(i)} 


be  endowed  with  the  inner  product 


<(^.4)(«.).0).(X.v(*.),S')>„  -  (i:-4>(i))(X->(.()L))  +  <D‘((>.D»^v>- 


where  the  SYmbol  D  is  used  here  and  below  to  denote  the  spatial 

d 

differentiation  operator  — .  The  space  V  together  with  the  inner 

dx 

product  < • , • >y  form  a  Hilbert  space  which  is  densely  and  compactly 
embedded  in  H. 

We  rewrite  the  system  (2.1)  -  (2.7)  as  the  abstract  second  order 
initial  value  problem  in  H 

(3.1)  Mu^^(t)  +  Cu^(t)  +  Ku(t)  =  F(t).  t  >  0 

(3.2)  5(u(t)  +  eu^(t))  -  0,  t  >  0 

(3.3)  u(0)  -  0  u^(0)  -  0 

in  the  states  u(t)  -  (y(t)  .u(t .«.)  ,u(t.  •  ))  .  The  operators  M  e  L(H)  , 
C:D  C  H  -►  H  and  K:D  C  H  -♦  H  are  given  by 

(3.4)  M(^,T|,4))  -  (mjj^  .m^Ti.pet))  , 

C(^,Ti,(t))  -  (Cjj(I;-ti),Cjj(ti-^)  -  Cj^I  D^()>0l),Cj5l  D*(t)). 

and 

K(;.Ti.(t))  -  (kjj(;-Ti).kg(n-^)  -  El  D’())(Jl),EI  D«(t)) 

where  D  -  {(^.n.^)  e  V:  <()  €  H‘(0.il)}.  For  each  t  >  0.  F(t)  - 
(0,f(t),0)  e  H,  6:  D  C  H  -♦  R  is  given  by  6((^,ti,(})))  -  D*(t)(!L)  and  e  - 

Cd/E. 

The  restrictions  C  and  K  of  the  operators  C  and  K  that  appear  in 
equation  (3.1)  above  to  Z/(5),  the  null  space  of  the  operator  6,  have 


natural  extensions  to  bounded  operators  from  V  (which  is  the  V-closure 
of  N(5))  into  V' ,  the  dual  of  V.  The  extensions  are  defined  in  terms 
of  the  bilinear  forms  c(-,-):  V  x  V  -•  R  and  k(  •  ,  •  )  :  V  x  V  -*  R  given  by 


(3.5)  (C(t))(s-)  -  c(J.S<)  -  Cjj( ?;-({> (Jl))(\-^^.(il))  +  Cj^I<D*(t),D*N^>Q 
and 

(3.6)  (K({>)(>*-)  “  k(J.H>)  =  kjj(;-ct)(iL))(X-rv(l))  +  El  <D*  0  .  >  q 

for  i  -  (^.<t)(il).<t))  €  V  and  ^  -  (X . ^^»(^.) , ^v)  6  V. 

The  finite  element  method  we  develop  below  could  be  derived  from 
standard  energy  considerations.  While  this  is  not  the  approach  we 
take,  it  is  worth  noting  that  the  usual  energy  expressions  can  be 
given  in  terms  of  the  forms,  operators  and  inner  products  defined 
above.  The  kinetic  energy  Is  given  by 

1 

^0  “  7  <Mu^(t).u^(t)>jj, 

the  potential  or  strain  energy  by 

1  ..  .. 

- k(u(t)  ,u(t)) 

^  2 

and  the  Rayleigh  dissipation  function  by 
1  .. 

J’o  “  “  c(u^(t)  ,u^(t))  . 

2 

Written  in  its  weak,  variational  or  distributional  form 

<Mu^^(t)  +  c(u^(t),(())  +  k(u(t).(t))  -  <F(t),<t)>g, 

t  >  0.  (p  €  V 


(3.7) 
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(3.8) 


u(0)  -  0 


u^(0) 


the  initial  value  problem  (3.1)  -  (3.2)  in  H  becomes  an  initial  value 
problem  in  V'.  If  we  assume  that  f  £  Li(0,T)  and  rewrite  (3.7),  (3.8) 
as  an  equivalent  first  order  vector  system ,  the  theory  of  abstract 
parabolic  systems  (see  [9],  [14])  yields  the  existence  of  a  unique 
mapping 

u  e  C([0,T];V)  n  HK(0.T):V)  n  C"([0.T];H)  n  H*((O.T);V') 

which  satisfies  (3.7),  (3.8).  If  we  are  willing  to  assume  further 
that  f  is  Holder  continuous  then  there  exists  a 

(3.9)  u  €  C([0,T];V)  n  G^((0,T];V)  n  CKtO,T];H)  n  C*((0,T];H) 

with  u(t)  +  eu^(t)  e  D,  t  >  0  which  uniquely  satisfies  the  initial 
value  problem  (3.1)  -  (3.3). 

In  order  to  demonstrate  the  convergence  of  the  approximation 
schemes  we  develop  below,  we  shall  require  a  somewhat  more  regular 
solution  to  the  initial  value  problem  (3.7),  (3.8)  than  either  of  the 
conditions  on  f  srated  above  can  guarantee.  In  addition  to  (3.9),  we 
shall  require  that  u  €  H*((0,T);V).  This  can  be  guaranteed  (see  [7]) 
if  we  assume  that  f  6  H‘(-t:,T)  for  some  x  >  0  with  f(t)  -  0,  t  <  0  and 
we  modify  our  original  mathematical  model  so  that 

(3.10)  F(t)  -  f(t)e,  t  e  [-x.T] 

for  some  0  -  (0,901). 9),  a  fixed  element  in  V.  We  note  that  with  9 


I 

I 


I 


■A 


$ 


» 

V" 

*•  H 

I 

•'“jI 

‘.‘-j 

v‘< 


.N 

•'V'ji 

I 

':j| 


-3 

•'J 

V.'. 

I 

i 


*  t  irf*  j.^  fM"  at  1*^ 


chosen  appropriately  in  V,  F  given  by  (3.10)  may  in  fact  represent  an 
improved  model  of  reality  when  compared  with  our  present  choice  of  F 
where  9  =  (0,1,0)  6  H. 

Central  to  our  approach  is  a  cubic  spline  based  Galerkin  approxi¬ 
mation  to  the  initial  value  problem  (3.7),  (3.3).  For  each  N  -1,2,- •• 

N  i  21 

let  A  denote  the  uniform  mesh  {0,  — ,  — ,  •••,  1}  on  [0,1]  and  let 

N  N 

N  N+1 

{B^}^  ^  denote  the  usual  cubic  B-splines  defined  with  respect  to  the 

N  N  2 

nodal  set  A  (see  [1-],  [13]).  Briefly,  each  B  is  a  C  function  on 

.] 

1  1 

[0,1]  which  is  a  cubic  polynomial  on  each  subinterval  [(k-1)— ,k— ], 

N  N 

N  ,  1  1,  ,  , 

k  =  l,2,-*-,N.  The  support  of  B  is  [ ( j-2)— ,  ( J+2)— ]  n  [0,l]  with 

j  N  N 

N  1^  N  1  N  1  N  1  N 

B  (V-)  -  4,  DB  (j-)  -  0,  B  ((J±l)-)  »  1,  and  DB  ((J±l)-)  =  +  -. 

JN  JN  JN  jNl 

N,N+1  N  N  N  N  N  N 

Defining  [0  }  by  0  =■  B  -2B  -2B  and  0  =■  B  ,  1-2 , 3 ,  •  •  •  ,N+1 , 

J  J-1  1  0  1  -1  J  i 

N  N 

we  „o.ve  0  (0)  -  D0  (0)  =  0,  J  -  1,2,---,N+1.  With  0  -  (1,0,0)  and 

j  J  0 


"N  ^  N  ^  N 

0  -  (0,0  (1),0  ).  J  -  1,2,---,N+1,  V 

«J  J  J 

dimensional  subspace  of  V. 


N  N+1 

span  10  }  is  an  N+2 

J  J“0 


The  Galerkin  equations  in  V  corresponding  to  (3.7),  (3.8)  for 
.  N 

u  (t)  e  V  are  given  by 


(3.11) 


(3.12) 


..N  .N  .N  .N  .N  .N  .N 

<Mu  (t),0  >  +  c(u  (t),0  )  +  k(u  (t),0  )  -  <F(t),0  >  , 

ttjHtJ  J  JH 


t  >  0,  J  -  0,1,2, • • • ,N+1 


-N 

u  (0)  -  0 


»N 

u  (0)  -  0. 
ti 


m 
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Setting 


.N  N+1  N 

^  (t)  -  w  (t)6  . 

J“0  J  J 


t  i  0, 


the  initial  value  problem  (3.11),  (3.12)  in  V  is  equivalent  to  the 
linear,  nonhomogeneous ,  second  order  N+2  -  vector  system 


d>w^  dw^ 

(3.13)  - (t)  +  - (t)  +  K^w^(t)  =■  F^(t), 

dt  *  dt 


t  >  0 


(3.14)  w  (0)  -  0 


-(0)  =  0 


N,  ^  N  N  T 

where  w  (t)  -  (w  (t),w  (t),’-*,w  (t))  .  The  entries  in  the  (N+2) 

0  1  N+1 

,  ^  N  N  N 

X  (N+2)  matrices  M  ,  C  and  K  are  given  by 


..N  .^N 

.N  .N 

c(0^,Bj), 


N  .N  ..N 


i,j  -  0 , 1 , 2 , • • • .N+1  respectively.  For  each  t  >  0  the  components  in 

N  N  .N  N 

the  N+2  -  vector  F  (t)  are  given  by  F  (t)  =  <F(t),0  >  -  f(t)0  (i) 

i  i  H  i 

or,  recalling  (3.10),  by 


F  (t)  -  f(t)<e.0  >  -  f(t)(e(j.)3  (i)  +  <e.0  >  ). 

i  i  H  i  i  0 


i  -  0,1,2, •• • ,N+1. 

We  consider  the  sequence  of  approximating  finite  dimensional  iden¬ 
tification  problems  which  consist  of  finding  6  Q  which  minimizes 

(3.15)  J-^Cq)  =  I - (t,iL:q)  -  z(t)  dt 

-’t,‘at*  I 

where  for  each  q  e  Q,  u^(t:q)  -  (y^(t;q)  ,u^(t  ,]l;q)  ,u^(t ,  •  ;q))  is  the 
unique  solution  to  the  initial  value  problem  (3.11),  (3.12)  in 
corresponding  to  q  -  (m,j,,E,CQ,mg,Cg,kjj)  6  Q.  In  actual  practice,  for 
a  given  q  €  Q,  J^(q)  is  computed  as 

J  (q)  -  +  4w\t;q)  +  wj^  (t:q)  -  z(t)|^dt 

-'to'  N-1  N  N-t-l  ' 

N  N  NT 

Where  w  (-jq)  =  (w  ( • ;q) , • • • ,w  (•;q))  is  the  unique  solution  to  the 

0  N-i-l 

N-i-2  -  vector  system  (3.13),  (3.14)  corresponding  to  q  €  Q. 

With  finite  dimensional  state  constraints,  the  solution  of  the 
N  estimation  problem  above  is,  at  least  in  principle,  routine.  For 
inverse  problems  which  are  closely  related  to  the  one  we  treat  here, 
our  earlier  numerical  studies  have  shown  that  satisfactory  results  can 
be  obtained  using  any  one  of  a  number  of  standard  computational 
techniques  for  least  squares  minimization  (for  example,  Newton's 
method,  conjugate  gradient,  steepest  descent,  Levenberg-Marquardt , 
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V’ 

V 

n' 


Our  fundamental  theoretical  result  is  that  each  of  the  approxi¬ 
mating  identification  problems  and  the  original  problem  have 
solutions.  Moreover,  we  show  that  the  solutions  to  the  approximating 
problems,  in  some  sense,  approximate  solutions  to  the  original 
problem.  We  require  the  following  lemma. 

T.ffTnTna.  3 . 1  Suppose  {q^}  C  Q  with  q*^  -*  q*^  as  N  -♦  m.  Let  u^(  •  ;q^) 
denote  the  unique  solution  to  the  initial  value  problem  (3.11),  (3.12) 
corresponding  to  q  and  let  u(-;q  )  denote  the  unique  solution  to  the 
initial  value  problem  (3.7),  (3.8)  corresponding  to  q*^.  If  u(-;q^)  € 
H*((0,T);V)  then 

2 

dt  -  0 
H 

as  N  -♦  00, 

Proof 

N 

For  each  N  -  1,2,*’-  let  P  denote  the  orthogonal  projection  of  H 
N 

onto  V  defined  with  respect  to  the  standard  inner  product  on  H, 
<',->jj.  Using  the  approximation  theoretic  properties  of  interpolatory 
splines,  it  is  not  difficult  to  show  that  (see  [3]) 

(3.17)  lim  |(P^  -  I)(2;,n.4))|  -  0 

N-»<»  I  I H 

for  each  (^,Ti,(t))  €  H  and  that 

(3.18)  lim  1(P  -  I)(t)|  -  0 

N-*oo  I  I V 


rT.  .N 


(3.16)  )  -  'i^^Ct;q°) 

^  0  UU  wU 


for  each  O  G  V. 


For  q  -  (m^.E.CQ.nijj.Cjj.kjj)  6  Q  it  is  immediatelY  clear  that  M, 
c(-,-)  and  the  operator  and  forms  defined  in  (3.4),  (3.5)  and 

(3.6)  respectively  depend  upon  q.  For  q*^  =  (m*^,E^,c®,ra^.c*^,k^)  €  Q 

T  D  H  H  H 

N  N  N  N  N  N  N 

and  q  -(m,E,c,m,c,k)€Qwe  adopt  the  shorthand  notation 
T  DHHH 

M°  -  M(q°).  C°(...)  -  C(q°)(.,.).  k°(..-)  -  k(q°)(..0.  -  M(q^) . 

c^(-,*)  -  c(q^)(*,*)  and  k^( • , • )  -  k(q^) (•,•)•  Similarly,  we  denote 
u(-:q  )  and  u  (•;q  )  by  u  and  u  respectively. 

From  (3.17),  the  assumption  that  u^  e  H*((0,T);V)  and  the 
inequality 

it  is  clear  that  we  need  only  to  consider  the  first  term  on  the  right 
hand  side  of  the  above  estimate. 

Letting  v^(t)  -  u^(t)  -  p\°(t)  for  t  ^  0,  (3.7),  (3.8),  (3.11), 
(3.12)  and  C  V  imply 


(3.19) 


;:n  :n.  ^  ;:Nn 

<h  ^tt’^  ’h  °  ^^t’^  )  +  k  (v  ,4)  ) 
-  <M^(I-P^)u°^,J^>2  +  <(M°  - 
.  c^((i-pN);o,;N) ,  o°(;o.;^)  - 


^  kN((i-pN);o.;N)  ^  t  >  o,  g 


21 


*&' 

$ 

K*;: 

•Wi 

I 

I 


'I***! 

*!v 

w 

.V. 


w 


I 

I 

s 

'•!i4 

**S4 
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*  t  ck°(i°.v»)  -  k^cio.v”))  -  k^cv^.v^). 

We  recall  that  Q  has  been  assumed  to  be  a  closed  and  bounded  subset  of 
R®  and  observe  therefore  that  the  forms  c®(-,-).  c^( ■ , • ) ,  and 

k^C • , * )  are  uniformly  bounded.  These  two  facts  together  with  the 
repeated  application  of  the  inequality 


<a.b>  i  |allbl  ^  Qt|al^  +  -  l^l^»  ^  *  0 

4a 


in  (3.21)  yield  the  estimate 


Jo  l-Li/-  * 


^  Jo  ^  “I’-C 


^  |v“| 


■  *  •  •  v^.-  V  ■.-  V  V 


where  Yq  is  a  positive  constant.  Choosing  a  >  0  sufficiently  small, 
we  find 

(3.22)  [  |v^  (s)|  ds  +  |v^(t)|  i  +  f  o, (s)ds  +  Y,  f  |v^(s)|  ds 

JqI  SS  Ijj  i  t  U  Jq  i  IJqI  S  ly 

where 

^2  >.2 
CTgCt)  -  Ygt  |(I-P^)’^°Ct)|^  +  |(I-P^)u°(t)|^ 


and  Yj_,  i  -  1.2,3  are  positive  constants  which  do  not  depend  on  N. 

Choosing  <|)^  -  v^(t)  €  in  (3.19),  arguments  similar  to  those 
used  above  (see  [2],  [3])  yield 


(3.23) 


I  I  ^ 

lim  v^(t)  =■  0 

N-o.  V 


for  each  t  6  [0,T].  Using  u^  €  H*((0,T);V),  (3.18)  and  an  application 
of  the  Gronwall  inequality  to  (3.22)  we  obtain  the  desired  result. 


Ve  note  that  we  also  obtain 
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(3.24) 


i*-"  -  ° 

N-*oo  V 


for  each  t  e  [O.T].  From  (3.23)  and  (3.24)  we  find  I u^(t : q^)-u(t ; q°) I 

'  V 

-♦  0  and  |u^(t  ;q^)-u^(t  ;q*^)  I  -  0  as  N  -►  oo  for  each  t  e  [O.T]. 

We  remark  that  it  is  the  L,  convergence  (more  precisely,  H 
convergence)  in  (3.16)  which  necessitates,  at  least  in  theory,  that  we 
be  provided  with  distributed  time  observations  (i.e.  observations 
which  are  continuous  in  time).  It  is  clear  from  (3.23)  and  (3.24) 
that  for  fits  based  upon  displacement,  velocity  or  slope,  time-sampled 
measurements  are  sufficient.  Of  course  when  the  approximating 
optimization  problems  are  solved,  the  integral  least  squares 
performance  indices  (3.15)  are  discretized.  Consequently,  in 
practice,  only  discrete  measurements  of  linear  acceleration  at  the  tip 
are  required. 


Theorem  5 . 1  Each  of  the  approximating  identification  problems  has  a 
“N  -N 

solution  q  .  The  sequence  {q  }  C  Q  admits  a  convergent  subsequence 
_N .  _N ,  _ 

{q  with  q'^-*q€(5asJ-*oo.  If  for  each  q  €  Q,  u(-:q),  the  unique 
solution  to  the  initial  value  problem  (3.7),  (3.8)  corresponding  to  q, 
is  an  element  in  H*((0,T):V)  then  q  is  a  solution  to  the  original 
Identification  problem.  In  addition,  the  limit  point  of  any  convergent 
subsequence  of  {q^}  is  a  solution  to  the  original  identification 
problem  as  well . 

Proof  Standard  continuous  dependence  results  for  linear  ordinary 
differential  equations,  the  fact  that  Q  has  been  assumed  to  be  a 
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I 


6  N 

closed  and  bounded  subset  of  R  and  the  form  of  J  are  sufficient  to 

conclude  that  a  solution  e  Q  to  the  approximating  identifica¬ 
tion  problem  exists.  Once  again  since  Q  is  a  closed  and  bounded  (and 

ft  — M 

therefore  compact)  subset  of  R  ,  the  sequence  {q  }  C  Q  admits  a 

_N .  _N  _N .  _ 

convergent  subsequence.  If  {q  ■J)  C  {q  }  with  q’^-*qtQasJ-»a3  and 
q  is  any  point  in  Q.  then  two  applications  of  Lemma  3.1  (the  second 
one  with  the  constant  sequence  {q})  yield 

N  _N  N 

J(q)  -  lira  J  ^(q  J)  i  lim  J  J(q)  -  J(q) 

J-*®  j-*® 


and  the  theorem  is  proved. 


Although  Theorem  3 . 1  above  guarantees  only  subsequential 
convergence,  in  all  test  and  simulation  examples  we  have  considered, 
we  in  fact  observe  the  convergence  of  the  sequence  {q^}  itself  to  the 
optimal  parameters  q.  Also,  it  is  not  difficult  to  verify  that  with 
only  minor  modification  (see  [2])  the  approximation  scheme  reported  on 
here  (together  with  the  convergence  theory  outlined  in  the  lemma  and 
theorem  above)  could  be  applied  to  inverse  problems  involving  the 
estmation  of  spatially  varying  parameters  (such  as  linear  mass  density 
p,  flexural  stiffness  El,  or  damping  coefficient  Cj^I)  which  appear  in 
the  equations  (2.1)  -  (2.4).  ¥e  note  of  course  that  when  either  El  or 
Cj^I  are  spatially  varying,  the  Euler-Bernoulli  equation  and 
corresponding  boundary  conditions  are  of  a  slightly  different  form 
than  those  given  in  (2.1)  -  (2.4)  (see  [3]). 


We  used  our  scheme  to  attempt  to  solve  the  inverse  problem  which 
was  posed  above  with  data  obtained  from  an  experiment  on  the  RPL 
structure.  We  report  on  our  findings  and  observations  here. 

All  computer  codes  were  written  in  Fortran  and  run  on  the  IBM 
3081  at  the  University  of  Southern  California.  The  approximating 
finite  dimensional  least-squares  minimization  problems  were  solved 
using  the  IMSL  implementation  of  the  Levenberg-Marquardt  algorithm 
(routine  ZXSSQ) ,  an  iterative  Newton's  method-steepest  descent  hybrid 
(see[2]).  The  second  order  N+2  -  vector  systems  (3.13),  (3.14)  were 
solved  (integrated)  in  each  iteration  (for  the  evaluation  of  and 
its  gradient)  using  Gear's  method  for  stiff  systems  (IMSL  routine 
DGEAR) .  The  integral  least  squares  performance  index  was  approximated 
by  a  discrete  sum  over  a  uniform  mesh  on  Ct.,tx].  The  integral  inner 
products  in  the  definitions  of  the  matrices  M^,  and  were 
computed  using  a  composite  two  point  Gauss-Legendre  quadrature  rule. 
The  second  time  derivative  of  w^,  or  generalized  acceleration, 

d^w^ 

— — ,  was  computed  using  a  second  order  centered  difference  on  the 
dt  * 

generalized  displacement, 

w^(t+A)  -  2w^(t)  +  w^(t-A) 

(4.1)  - (t)  - - 5 - . 

dt  *  A*^ 

We  found  this  to  be  a  somewhat  more  stable  method  for  computing 
acceleration  (an  unbounded  measurement)  than  was  a  first  order 
centered  difference  on  the  generalized  velocity. 
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(4.2) 


N  N 

dw”  A  dw  ,  A 

- (t  +  -) - (t  -  -) 

dt  2  dt  2 

A 


Either  of  the  time  differencing  formulas  (4.1)  or  (4.2)  proved  to  be 
significantlY  more  stable  than  using  the  differential  equation  (3.13) 

N 

directly  to  compute  ~ ^ -( t )  via  an  inversion  of  M  .  As  to  why  this 

was  so,  we  can  only  offer  the  conjecture  that  the  time  differencing 

provided,  at  least  to  a  certain  extent,  some  filtration  of  the  signal. 

Before  turning  our  attention  to  the  experimental  data,  we  tested 

our  scheme  with  simulated  data.  "True"  values  for  the  unknown 

parameters  c^  (actually  CqI),  and  k^  were  chosen  and  a  quintic 

spline-based  semi-discrete  Galerkin  scheme  applied  to  the  initial 

value  problem  (3.7),  (3.8)  was  used  to  generate  data. 

Setting  p  .03,  m  ,J,-  .15,  El  -  80.0,  1-4.0  and 

f(t)  -  {  0  i  t  i  0.05 

0.0  0.05<tS5.0, 

the  fit  was  carried  out  based  upon  observations  of  linear  acceleration 
at  the  tip  at  times  t^  -  .11,  i  -  2. 3, •••,50.  We  note  that  this  is 
equivalent  to  taking  t,  -  .1,  ti  -  5.0  and  using  a  standard  rectangle 
rule  with  uniform  mesh  spacing  .1  to  discretize  the  integral  appearing 
in  the  definition  of  the  least  squares  performance  index  J^.  The 
initial  estimates  c^I  -  .0035,  m^  -  .035 
start  the  iterative  optimization  procedure.  In  (4.1),  A  was  taken  to 
be  .1.  Our  results  are  summarized  in  Table  4.1  below. 


,  and  kjj  ”  were  used  to 


5  iL.' 


2 

.037537 

.039471 

.003428 

.298626 

2.57x10“^ 

3 

.066997 

.039485 

1  .003907 

.298875 

4.37x10"^ 

4 

.005063 

.039777 

1 

.003997 

. 299455 

5.06x10"^ 

5 

.005667 

.039899 

.003971 

.299787 

7.66x10"^ 

6 

.005049 

.040035 

.004006 

.300087 

4.63x10“® 

True 

value 

.005000 

.040000 

.004000 

.300000 

Initial 

Estimate 

.003500 

1 

1 

.035000 

.003500 

.400000 

Table  4.1 


The  experiment  which  we  describe  below  was  carried  out  for  us  on 
the  RPL  structure  by  Dr.  Michel  A.  Floyd,  formerly  of  the  Control  and 
Flight  Dynamics  Division  of  the  Charles  Stark  Draper  Laboratory  and 
the  Department  of  Aeronautics  and  Astronautics,  MIT. 

The  air  bearing  table  was  clamped  so  that  the  central  hub  could 
not  rotate.  The  thruster  lines  for  one  of  the  active  appendages  was 
set  to  300  psi  and  the  thruster  was  fired  for  .05  seconds  (50  milli¬ 
seconds).  With  the  appendage  initially  at  rest,  the  firing  of  the 
thruster  was  assumed  to  have  begun  at  time  t  -  0.  Linear  acceleration 
at  the  tip  was  observed  over  the  time  interval  0  to  5  seconds.  With  a 
sampling  period  of  .005  seconds  (5  milliseconds)  a  total  of  1000 
measurements  were  recorded.  The  data  is  plotted  in  Figure  4.1  below. 
The  scale  factor  for  the  accelerometer  is  5  volts/g  (g  -  32  ft/sec*)- 
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Figure  4.1 

The  noticeably  higher  frequency  (-  14  Hz)  component  of  the  data 
is  a  torsional  mode  of  the  arm  excited  by  the  motion  of  the  thruster 
valve  mechanisms  and  Inertial  and  elastic  forces  applied  to  the  tip  of 
the  arm  by  the  nitrogen  supply  hose.  The  opening  or  closing  of  the 
solenoidal  valve  in  the  thruster  generates  an  inertial  force  which 
acts  as  a  torque  on  the  tip  of  the  arm.  Consequently,  torsional  modes 
are  excited.  Also,  in  addition  to  modifying  transverse  bending 
characteristics,  since  the  hose  is  attached  to  the  top  of  the  arm,  its 
horizontal  motion  will  tend  to  generate  torques  which  have  a 
"twisting"  effect.  Although  the  accelerometer  is  mounted  at  the 


center  of  the  arm  (and  therefore  on  a  nodal  line  of  the  longitudinal 


torsional  modes,  if  we  assume  vertical  symraetrY) .  as  the  arm  twists, 


the  accelerometer  picks  up  a  component  of  the  earth's  gravitational 


force.  Since  the  first  torsional  mode  has  a  much  higher  frequency 


than  either  of  the  first  two  flexible  modes  (.75  Hz  and  7.5  Hz,  as 


identified  from  an  FFT  of  the  data)  and  since  it  is  rapidly  damped,  we 


neglected  its  contribution  to  the  accelerometer  signal,  treating  it  as 


white  noise,  and  left  it  unmodeled.  A  detailed  discussion  of  the 


causes  of  the  excitation  of  the  torsional  modes  and  its  effect  on  the 


transverse  bending  characteristics  of  the  active  appendages  can  be 


found  in  [6] . 


The  physical  characteristics  of  the  structure  are  as  follows. 


The  arm  is  made  of  aluminum  and  is  4  feet  in  length,  6  inches  in  width 


and  .125  inches  in  height.  From  this  we  obtain  X.  »  4.0  ft,  p  =  .027 


slug/ft  and  I  =  4.71  x  10~®  (ft)"^.  The  theoretically  predicted  value 


for  E  is  15.84  X  10®  lb/(ft)^.  The  mass  of  the  thruster  assembly  was 


determined  to  be  m,_  -  .149  slug.  From  the  calibration  table  in  [6], 


we  find  that  a  hose  pressure  of  300  psi  is  equivalent  to  a  force  of 


297  lb.  We  set  therefore 


f(t)  -  { 


0.297  lb 


0  S  t  5  0.05 


0.05  <  t  ^  5.0 


To  serve  as  a  basis  for  comparison,  we  neglected  the  hose  effects 


and  structural  damping  (i.e.  we  chose  c^^  -  m^  -  c^  -  k^j  -  0)  and  used 
the  standard  Euler-Bernoulli  model  with  the  parameters  p,  E,  I  and  m,^ 
and  input  f  as  specified  above  to  generate  the  plot  of  linear 


acceleration  at  the  tip  given  in  Figure  4.2. 
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Figure  4.2 


The  plot  was  obtained  by  Integrating  the  initial  value  problem  (3.13), 
(3.14)  with  N  -  4  and  then  using  (4.1)  to  compute  the  acceleration  at 

3  *  u  3  *  u^ 

the  free  end.  The  residuals  ( — ^(t.X.) - (t.i))  over  the  time 

3ti  *  3 1  * 

interval  [0,5]  are  plotted  in  Figure  4.3.  The  sum  of  the  squares  of 
the  residuals  (at  intervals  of  .1  seconds)  was  found  to  be  3.03. 
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Figure  4.3 


Using  the  data  on  the  interval  3.0  to  5.0  (where  the  contribution 
from  the  torsional  inodes  has  been  significantly  damped)  with  a 
sampling  period  of  .1  seconds  we  used  our  scheme  with  N  -  4  to  obtain 
optimal  estimates  for  the  coefficient  of  viscosity  Cj^  and  the  hose 
parameters  m^,  c^  and  kjj.  In  the  set  of  runs  we  are  about  to  describe 
the  values  of  E  and  m^,  were  held  fixed  at  their  theoretically 
predicted  values.  A  rough  calculation  based  upon  "matching"  the  first 
two  observed  natural  frecjuencies  of  the  data  with  the  first  two  modal 
frequencies  of  the  model  was  used  to  obtain  a  crude  initial  estimate 


for  the  ratio  kjj/nijj.  Then,  using  our  scheme  to  minimize  over  the 
parameters  mj^  and  only .  we  obtained  the  optimal  values  shown  in 
Table  4.2  below.  Integrating  the  system  (3.13),  (3.14)  over  the  time 
interval  [0,5]  with  m„  and  k„  set  to  the  values  in  the  table  and  c„  - 

n  il  u 

=  0  the  sum  of  the  squares  of  the  residuals  (at  intervals  of  .  1 
seconds)  was  found  to  be  .73. 


mjj(slug) 

kjjdb/ft) 

.039269 

.339935 

Table  4.2 


Next,  holding  m^j  and  kjj  fixed  at  the  values  shown  in  Table  4.2,  a 
search  on  c^  was  carried  out  (the  initial  estimate  for  c^  was  taken  to 
be  zero  and  c^  was  held  fixed  at  zero).  Then  using  the  resulting 
values  of  m^^,  c^  and  k^  as  initial  estimates,  a  fit  over  all  three 
parameters  was  performed.  The  result  is  shown  in  Table  4.3.  The  sum 
of  the  squares  of  the  residuals  was  found  to  be  .728. 


mjj(slug) 

Cjjdb  •  sec/ft) 

kjjdb/ft) 

.043431 

.004056 

.351385 

Table  4.3 


Continuing  to  use  the  same  procedure  to  generate  "start  up"  values,  we 


eventually  used  our  scheme  to  search  over  all  four  parameters  o^,  m^^ 

c„  and  k„  simultaneously  obtaining  the  values  given  in  Table  4.4  and 
11 

the  fit  plotted  in  Figure  4.4.  The  residuals  are  plotted  in  Figure 
4.5.  The  sum  of  their  squares  was  computed  to  be  .70. 


c^Clb-sec/Cft)^) 

mjj(slug) 

c„(lb • sec/ft) 

ll 

k„(lb/ft) 

il 

127.40 

.0801 

.007804 

.412977 

Table  4.4 
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Figure  4.4 
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Figure  4.5 


In  designing  a  controller  for  the  RPL  experiment,  Floyd  in  [6] 
used  model  adjustment  to  tune  a  simple,  undamped,  cantilevered  beam 
with  tip  mass  model  for  the  active  arms  (i.e.  the  arms  with  the  hoses) 
of  the  structure.  He  used  the  following  procedure.  The  air  bearing 
table  was  locked  in  a  stationary  position.  With  the  hose  depressur¬ 
ized,  an  impulsive  force  was  applied  to  the  beam  and  linear  accelera¬ 
tion  at  the  tip  was  measured  and  recorded.  Based  upon  the  physical 
assumption  that  with  the  hose  depressurized,  the  presence  of  the  hose 
serves  only  to  add  mass  to  the  tip  of  the  arm,  the  parameter  m,j,  was 


A.-'-. 


adjusted  so  that  the  first  mode  or  frequency  of  the  model  agreed  with 
the  first  observed  cantilever  mode  (obtained  via  an  FFT)  of  the  data. 
Then,  with  the  hose  pressurized,  the  same  experimental  procedure  was 
carried  out.  This  time  however,  the  modulus  of  elasticity  E  of  the 
beam  was  adjusted  to  compensate  for  the  variation  in  stiffness  which 
results  from  the  presence  of  the  hose.  The  adjusted  values  of  the  tip 
mass,  m,j,,  and  modulus  of  elasticity,  1.  obtained  by  Floyd  are  given  in 
Table  4.5  below. 


m,j,  (slug) 

E  (lb/(ft)^) 

.254 

17.31  X  10® 

Table  4.5 


We  integrated  the  system  (3.13),  (3.14)  using  the  adjusted  values  of 
m,j,  and  E  given  in  the  table  (and  =■  =•  Cjj  =  =  0)  and  obtained 

the  plot  shown  in  Figure  4.6.  The  corresponding  residuals  are 
plotted  in  Figure  4.7.  The  siira  of  the  squares  of  the  residuals  was 
computed  to  be  5.1. 

Starting  with  the  same  basic  model,  we  used  our  scheme  to 
determine  the  values  of  m,j,  and  E  which  minimize  the  sum  of  the  squares 
of  the  residuals  over  the  time  interval  [3.0,  5.0]  with  a  sampling 
period  of  .1  seconds.  Taking  the  theoretically  predicted  values  of  m,^, 
and  E  (m,j,  -  .149  slug,  E  -  15.84  x  10®  lb/(ft)^)  as  start  up  values 
for  the  optimization  routine  yielded  the  results  given  in  Table  4.6. 


I 
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The  corresponding  fit  and  residuals  are  plotted  in  Figures  4.8  and  4.9 


Figure  4.7 
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Figure  4.8 
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Figure  4.9 

In  summary ,  we  have  seen  that  analysis  of  the  RPL  experimental 
data  can  be  carried  out  in  several  ways  with  a  number  of  different 
models.  Our  techniques  can  be  used  to  provide  reasonable  fits  of  the 
data  to  models  with  or  without  hose  and/or  beam  damping.  Even  if  one 
attempts  to  leave  the  physics  of  the  hose  -  beam  dynamic  interaction 
unmodeled  and  perform  "model  adjustment"  (by  adjusting  the  values  of 
the  tip  mass  m^  and  beam  modulus  of  elasticity  E),  our  estimation 
techniques  provide  a  much  better  fit  than  that  obtained  using  "modal 
matching"  methods  common  in  engineering  practice. 

One  of  the  primary  objectives  of  our  effort  here  was  to 
demonstrate  the  efficacy  of  our  scheme  and  in  particular,  to  assess 
its  effectiveness  when  provided  with  actual  experimental  data.  While 
we  are  pleased  with  the  results  obtained  for  the  RPL  data,  we  are 
careful  to  point  out  that  to  provide  a  fair  and  complete  evaluation  of 
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the  usefulness  of  our  models  for  the  RPL  experimental  structure,  a 


more  complete  and  in-depth  study  involving  extensive  experimental  work 


and  statistical  analysis  would  necessarily  be  required. 
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